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Abstract 

Transport of low energy neutrons associated with the galactic cosmic ray 
cascade is analyzed in this dissertation. A benchmark quality analytical algorithm 
is demonstrated for use with BRYNTRN, a computer program written by the High 
Energy Physics Division of Nasa Langley Research Center, which is used to design 
and analyze shielding against the radiation created by the cascade. BRYNTRN uses 
numerical methods to solve the integral transport equations for baryons with the 
straight-ahead approximation, and numerical and empirical methods to generate 
the interaction probabilities. The straight-ahead approximation is adequate for 
charged particles, but not for neutrons. As Nasa Langley improves BRYNTRN to 
include low energy neutrons, a benchmark quality solution is needed for 
comparison. The neutron transport algorithm demonstrated in this dissertation 
uses the closed-form Green s function solution to the galactic cosmic ray cascade 
transport equations to generate a source of neutrons. A basis function expansion 
for finite heterogeneous and semi-infinite homogeneous slabs with multiple energy 
groups and isotropic scattering is used to generate neutron fluxes resulting from 
the cascade. This method, called the Fn method, is used to solve the neutral 
particle linear Boltzmann transport equation. As a demonstration of the algorithm 
coded in the programs MGSLAB and MGSEMI, neutron and ion fluxes are shown 
for a beam of fluorine ions at 1000 MeV per nucleon incident on semi-infinite and 
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finite aluminum slabs. Also, to demonstrate that the shielding effectiveness against 
the radiation from the galactic cosmic ray cascade is not directly proportional to 
shield thickness, a graph of transmitted total neutron scalar flux versus slab 
thickness is shown. A simple model based on the nuclear liquid drop assumption is 
used to generate cross sections for the galactic cosmic ray cascade. The 
Endf/B V database is used to generate the total and scattering cross sections for 
neutrons in aluminum. As an external verification, the results from MGSLAB and 
MGSEMI were compared to Anisn/Pc, a routinely used neutron transport code, 
showing excellent agreement. In an application to an aluminum shield, the F N 
method seems to generate reasonable results. 
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Chapter 1 
Introduction 

The advance of man and machines into space presents concerns about 
damage caused by radiation. This damage can disrupt electronic and biological 
systems to the point of failure. To protect against radiation damage, three 
mechanisms are available: distance, time, and shielding. Providing distance 
between the radiation source and the systems in a space environment is expensive 
because of size and weight limitations. Time considerations are not applicable 
because the spacecraft is exposed continuously to radiation sources. Therefore 
shielding must be used as the primary protection mechanism. 

The radiation of concern in the space environment is highly energetic 
(billions of electron volts per nucleon) heavy ions and their secondary radiation. 
The primary ions originate from deep space sources or from our sun and are called 
galactic cosmic rays. The secondary radiations are fragments resulting from 
collisions of galactic cosmic rays with the material used to shield the systems in 
question. As these fragments collide with shield material, smaller fragments are 
created. Eventually, the fragments generated are neutrons and protons. The 
resultant radiation shower created by this phenomenon is called the galactic 
cosmic ray cascade. The prediction of the radiation dose resulting from this source 
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is very important for the proper design of shields for occupants and electronic 
components. 

The nuclear power industry uses mass to shield systems. This mass 
attenuates particles by multiple scattering collisions or by absorption. The 
particles being shielded usually do not fragment and create a secondary radiation 
shower. Using current nuclear power industry shield designs, large quantities of 
mass would be placed into orbit which is expensive and time consuming. To 
reduce mass in a spacecraft, the shielding used for galactic cosmic rays must 
perform other tasks, such as being paxt of a pressure vessel, a micro-meteor shield, 
or a structural component. A spacecraft’s less massive shielding must not only 
shield the internal systems from galactic cosmic rays, it must perform its other 
design tasks without creating a secondary radiation field that is worse than the 
original radiation field it was meant to shield. 

Through various efforts, personnel at the High Energy Physics division of 
the NASA Langley Research Center have developed a computational algorithm 
which employs numerical methods and various empirical and analytical nuclear 
physics computations to predict the radiation dose from the galactic cosmic ray 
cascade. The transport model used in this algorithm, while being realistic for 
charged particles, is not realistic for low energy neutrons. The work presented in 
this dissertation addresses the neutron problem by coupling an analytical solution 
to the galactic cosmic ray cascade with an algorithm that treats low energy 
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neutrons realistically. 


1.1 Motivation 

Various solution methodologies axe available to predict the radiation dose 
from the galactic cosmic ray cascade. The computer program that Nasa has 
created is called BRYNTRN (Reference [1]). It is a compromise between 
computational accuracy and ease of use in design tasks. The use of the Monte 
Carlo method or perturbation theory to design shielding can be highly accurate at 
the expense of long turn-around times and extensive computational resources. An 
Sn based deterministic method utilizing SMART scattering (References [2] and [3]) 
could be used to solve the problem in one, two, or three dimensions. While being 
faster than the Monte Carlo method and almost as accurate, it still consumes large 
amounts of computer resources. Resource usage could be reduced by limiting the 
number of particles treated and the number of dimensions. Of course, a reduction 
in utility and accuracy is also incurred. The numerical solution method employed 
by BRYNTRN using the integral transport equation requires a small fraction of the 
computational resources needed by Monte Carlo and is comparable to the S N 
method. Unfortunately, the physics approximations for neutrons reduce the 
applicability of the results in comparison to the other methods. In addition, 
verification of the BRYNTRN program, using a quality benchmark, is required to 
determine the accuracy of the dose calculation. 

The BRYNTRN program numerically solves a set of coupled integral 
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transport equations for the initial ion distribution and each subsequent fragment 
distribution created (including protons and neutrons) in energy and space to 
determine dose rates. The interaction and particle generation probabilities are 
determined from analytical and empirical nuclear physics calculations. The results 
for relatively low energy (100 to 400 MeV) neutrons and protons on tissue agrees 
well with various studies using three dimensional Monte Carlo codes 
(References [4] and [5]), but with more efficient use of computer resources. 

To verify the numerical solution technique used in BRYNTRN, a closed-form 
analytical solution to the galactic cosmic ray cascade equations wets developed in 
Reference [6] and used as a benchmark. The two methods generated virtually 
identical results over a wide range of input parameters. This verified the solution 
technique, but not the physics models on which they were based. To continue the 
development process, the models used in BRYNTRN must be expanded to 
encompass more realistic physics. To verify these new models, new benchmarks 
must be developed which is one purpose of the effort. 

The current BRYNTRN model treats neutrons as monoenergetic and 
monodirectional; therefore, they are not included in the overall dose calculation. 

In this dissertation, the closed-form analytical solution to the galactic cosmic ray 
cascade is coupled to an analytical neutral particle transport solver to generate a 
benchmark quality analytical solution to the cascade that includes a multiple 
energy group (multigroup) and angular dependent neutron description. With this 
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new algorithm as a benchmark, a better neutron model can be introduced into 
Bryntrn, tested, and verified. This will allow neutrons to more easily be included 
in the overall dose calculation. 


1.2 Background 

Radiation is very damaging to the systems of a spacecraft. Radiation 
damage is the destruction of the structure of system components by energetic 
particles. An extensive effort has been expended in formulating techniques to 
predict the damage caused by radiation, called dose, in various situations 
encountered in space. These dose predictions allow the design and building of 
shields to protect systems that axe sensitive to radiation. 

Analysis techniques for solving shielding problems are widespread in the 
nuclear industry; however, only neutrons of no more than 20 MeV, alpha particles, 
protons, beta particles, and gamma rays have been studied in great detail. Fission 
fragments in nuclear reactors, which axe heavier and more energetic than the above 
particles, are assumed to be confined to the nuclear fuel under normal operating 
scenarios. Fission fragment behavior is closer to the situation addressed in this 
dissertation than the other particles mentioned above, but the analysis techniques 
are almost entirely empirical. The main issue regarding fission fragment behavior 
centers on how long a fuel pellet can be exposed in a reactor before damage is 
extensive enough to warrant concern. The industry estimates the problem of 
fission fragment damage to the fuel pellet by using the integrated exposure, or 
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burnup, of the fuel pellet as a basis to measure other performance parameters. 
This large body of knowledge is of little use in generating a benchmark quality 
solution that contains the necessary detail for the galactic cosmic ray cascade. 

Since prior knowledge about charged particle interaction with matter 
cannot be used and built upon, a first principles approach is required to construct 
models that can be used to generate the necessary shield design. Experience has 
shown that the linear Boltzmann particle transport equation (hereafter, referred to 
as the Boltzmann equation) can be used as an approximate mathematical model 
for most problems concerning particle interactions. The full Boltzmann equation is 
difficult, if not impossible, to solve analytically. To create a tractable 
mathematical problem, various physical approximations and restrictions are 
required to simplify the full Boltzmann formulation. Previous work on the galactic 
cosmic ray cascade has shown that the Boltzmann equation without angular 
deflection is appropriate for high energy, heavy ion and subsequent proton fluxes 
(References [7], [8], and [9]). 


1.3 Specific Objectives 

As indicated above, the current models used in BRYNTRN should be 
expanded to include energy and angular dependent neutron distributions because 
the straight-ahead approximation, no angular deflection, does not accurately 
describe neutron motion in general. This is accomplished by coupling an existing 
analytical solution to the galactic cosmic ray cascade to an analytical neutral 
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particle transport solution. Two geometries are utilized in the neutral paxticle 
transport method: homogeneous semi-infinite and heterogeneous finite slabs. Both 
geometries are used to solve the galactic cosmic ray cascade problem. 

Homogeneous finite slabs can be used because the straight-ahead approximation 
does not allow the back-propagation of information about a boundary. The 
particles incident on a slab travel straight into the slab, but are not allowed to 
scatter, so they never travel back towards the slab boundary where they entered. 
Therefore, the particles never detect a boundary until they cross it, then they 
never re-cross that boundary. 

The closed-form analytical solution to the galactic cosmic ray cascade is 
determined by applying the Laplace transform to the appropriate Boltzmann 
equations. A set of ordinary, first order differential equations is generated and 
solved. Using a partial fraction expansion for the resultant solution, the inverse 
Laplace transform can be performed analytically. This determines the number, 
position, and energy of all the ions in transport media. To determine the number, 
position, and energy of the neutrons created, the Boltzmann equation is solved for 
neutrons without considering their motion in space and energy. The result is used 
in the source for the neutral particle algorithm. 

The neutral particle Boltzmann transport equation for one-dimensional 
heterogeneous finite and homogeneous semi-infinite media with multiple energy 
groups and isotropic down scatter is used as the model for the neutrons resulting 
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from the galactic cosmic ray cascade. This formulation contains an isotropic 
distributed neutral particle source. Since the neutron source described above is 
monodirectional, the collision operator from the Boltzmann equation is applied to 
the monodirectional source to create an isotropic source. The Boltzmann equation 
is transformed into two singular integral equations obtained by essentially taking 
the Laplace transform with respect to both the positive and negative directions 
and restricting the complex transform variable to be on the cut along the real axis. 
From the definition of a principle value of a Cauchy type integral, a set of 
Fredholm integral equations and constraints follow. These equations are solved 
using the standard method of expanding the solution in basis functions 
(Reference [10]) originally developed by C. E. Siewert in References [11] 
through [14]. To produce an answer relatively free from truncation error, a post 
processor based on integral transport theory provides the ultimate angular fluxes. 

The neutral particle transport algorithm can be decoupled from the galactic 
cosmic ray cascade source to provide a general benchmark quality neutral particle 
solution algorithm. A beam or isotropic source incident on the left face of the slab 
or a distributed source is available for heterogeneous finite or homogeneous 
semi-infinite slabs. 

The results generated by the neutral particle transport programs are 
angular and scalar fluxes. When the programs are in the galactic cosmic ray 
cascade source mode, ion fluxes and neutron source values are generated. This 



26 

information is presented in printable data files and data files that can be used to 
generate two-dimensional and three-dimensional plots. 

In Chapter 2, the physical theory and the associated mathematical models 
axe described. Once the mathematical models have been developed, the numerical 
methods used to solve the model equations are described in Chapter 3. Since the 
neutral paxticle transport program has been created for this dissertation, it is 
verified as described in Chapter 4. Once the program is verified, it is applied to 
the galactic cosmic ray cascade, and the results are discussed in Chapter 5. 
Concluding remarks and observations are contained in Chapter 6. Details of the 
derivation for the analytical neutron transport model can be found in Appendix A. 
A user’s manual for the programs is found in Appendix B. 
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Chapter 2 
Theory 

The galactic cosmic ray cascade is a complicated particle interaction that 
can be treated using classical physics. This chapter outlines physical and 
mathematical arguments required to simplify the problem so it can be solved using 
analytical techniques. To perform this task, a physical model is described. This 
physical model is then transformed into a mathematical model. 

The physical model details the interactions between the incident ion and 
the target material. Because of its speed, the ion is treated as a randomly 
configured mass of neutrons and protons traveling in a straight line. From the ion 
point of view, the target material is an evenly distributed proton and electron 
mass with distinct interaction centers. The evenly distributed mass acts 
continuously to slow the ion down though electrostatic interactions. Eventually, 
the ion slows and stops or encounters an interaction center and fragments creating 
smaller ions. The fragmenting continues until protons and neutrons are created. 
The protons act like the ions and eventually slow down and stop. The neutrons 
interact with the target material through scattering and absorption. 

The mathematical model used to represent the physical model is based on 
classical statistical mechanics through the linear Boltzmann particle transport 
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equation. This imposes physical and mathematical requirements on the solution; 
however, the assumptions required to use this equation generally exclude quantum 
mechanical and single particle effects through using a statistically large number of 
classical point particles. These assumptions do not effect the physical model 
enough to warrant concern. 

This chapter details the physical and mathematical models plus solution 
methods for the mathematical model equations. Once a solution method exists, 
numerical techniques can be used to generate numerical values which must then be 
verified. 


2.1 Physical Models 

The physical models required to describe this problem in a detailed manner 
do not represent the physics of the entire problem, but describe the physics 
sufficiently well to generate useful results. Definitions of various entities are needed 
to establish a basis for building a physical model. There are standard assumptions 
associated with the statistical, non-relativistic Boltzmann equation solutions which 
are described in detail and justified in Reference [15], Section 1.4. The first 
assumption is that all particles are treated as point particles; that is, a particle is 
completely described by its position and momentum. Spatial and temporal scales 
of the particles are sufficiently large to exclude quantum mechanical effects. Large 
numbers of each particle are present so deviations from the expectation value of 
the number densities can be ignored. The transport medium (target) does not 
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change on time scales of importance to particle transport. The density of particles 
relative to the density of the medium is small so that particle-particle interactions 
can be neglected. 

2.1.1 Galactic Cosmic Ray Interaction 

This model represents the interaction of heavy, high energy ions (galactic 
cosmic rays) with a shield or target material. The two modes of interaction 
considered are ions having coulombic collisions with target nuclei and slowing 
down, and ions colliding by direct impact with target nuclei creating nuclear 
fragments. Other modes of interaction exist and are occurring in galactic cosmic 
ray interactions with nuclei, but these are less significant and ignored in this 
treatment to reduce the complexity of the mathematical model. Detailed 
definitions of these two interactions are needed prior to generating mathematical 
models. 

As a heavy fast ion travels through matter, it loses energy by interacting 
with the matter through electronic excitation. This excitation is caused by the 
protons of the ion interacting with the protons and electrons of the target nuclei. 

In this process, energy and momentum are transferred to the nuclei from the ion 
by the electrostatic forces between the ion and the nuclei. This results in the ion 
slowing down and losing energy with each interaction. For this model, the gain in 
energy by the target nuclei is ignored. Since the target nuclei surround the ion, the 
media is considered continuous, so the slowing down process is considered 
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continuous. The ion is traveling at great speeds, so its direction is not changed 
substantially by the small pushes and pulls of the surrounding target nuclei. 
Therefore, the ion path is assumed to be a straight line. 

If the ion trajectory intersects a target nucleus, then a direct impact 
collision occurs called an abrasive collision. The result of an abrasive collision is 
fragmentation of the ion into other nuclear particles in a statistically random 
manner while conserving charge, energy, and momentum. To reduce the 
complexity of the model, the fragmentation of the target nuclei is ignored. 

These definitions can be restated by adding physical assumptions to the 
original assumptions required for the Boltzmann equation. Because the target 
material protons do not change the direction of the ions and fragments, the 
straight-ahead approximation is used to model the particle motion. Energy is lost 
by the ions in a continuous manner so the continuous slowing down approximation 
is used. All interactions axe statistically random while conserving charge, energy, 
and momentum. Fragmentation and energy gain of the target nucleus is neglected. 
From these assumptions and approximations, an appropriate mathematical 
interaction model can be formed. 

2.1.2 Neutral Particle Transport 

This model represents the interactions of neutrons with the target material. 
Since there is no coulombic charge associated with the neutron, only direct 
interactions are modeled. These interactions are either scattering collisions or 
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absorptions. The probability of each depends on the target material and the 
energy of the neutron. The Fn method, a robust mathematical and computational 
solution technique, is chosen to solve the neutron- Boltzmann transport equation 
for one-dimensional heterogeneous finite and homogeneous semi-infinite media 
with multiple energy groups and isotropic down scatter. An isotropic source is 
used to couple this model to the galactic cosmic ray cascade model. 

2.2 Mathematical Formulation 

The previous section described the physical models that require translation 
into a mathematical formulation. From the physical assumptions outlined, a 
statistical mechanics approach is taken. Therefore, a particular form of the 
Boltzmann equation is used for each model. This equation was originally used to 
describe the behavior of dilute gases; however, with appropriate modifications, it 
can be used to describe the behavior of dilute charged and neutral particles 
interacting with an appropriate medium. Appendix A describes the derivation of 
the Boltzmann equation for the Fn method. 

2.2.1 Boltzmann Equation 

There are various ways to derive the Boltzmann equation. The most 
fundamental is a physically based heuristic method (References [16] and [15]). This 
involves describing the various physical interactions and movements of particles in 
an infinitesimal phase volume. The Boltzmann equation can also be derived 
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directly from the statistical mechanical Liouville Equation (Reference [17]). The 
result of these derivations including only scattering and absorption interactions is 

'1 d 1 

• V r + a-V v + <r(r, E) $(r, 17, E, t) = 

= /°°dE' [ dO' <r(r, E') f(r;f2',E' -► fi, E) $(r, tl',E ’,t) + ( 2 - 1 ) 

JO J Air 

+ Q(r, E, t), 

where, 

v - Particle velocity, 
v - Particle speed = |v|, 

a - External force term acting on the particle = where, m is the 
particle mass, 

r - Position vector (three components), 

- Direction unit vector (two components) from the particle momentum 
= where, p is the particle momentum, 

E - Particle kinetic energy = Imu 2 , 

t - Time, 

cr(r, E) - Total interaction probability at position r and energy E per 
particle length of travel, 

- Particle angular flux [u7V(r, fi, E, t)] which is the speed times 
the particle density in particles per phase volume, 
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f(r; ft', E' — ► ft, E)dftdE - Probability density for particle transfer from ft' 
and E' to ft and E in dft and dE, 

Q(r, ft, E, t) - External particle sources. 

To determine a unique solution to this equation, initial and boundary conditions 
are required. The existence and uniqueness of a solution to the Boltzmann 
equation are not discussed here, but Reference [18] gives various proofs for the 
existence and uniqueness of a solution. A general, analytical solution to 
equation (2.1) at this point is impossible to obtain; however, the physical model is 
restrictive in various respects to allow simplifications to be made to the Boltzmann 
equation allowing a solution. A balance is struck between the physical restrictions 
and the mathematical simplifications in order to maintain a worthwhile problem 
relating to spacecraft shielding. 


2.2.2 Galactic Cosmic Rays 

As already indicated, the Boltzmann equation will be adapted to describe 
the galactic ray cascade by applying assumptions and approximations to reduce its 
complexity. If steady state conditions, straight-ahead motion, continuous slowing 
down, and multiple species considerations are assumed, then equation (2.1) 
reduces to the following set of equations for each ion species of charge number j 


d_ 

dr 


5e S j ( e ) + c ‘ 


0,(r,E) 


J 


z 

*=>+! 


M hk 


&k 0fc(T, E), 


(2.2a) 
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with boundary conditions of 

^(O.E) = /i(E), (2.2b) 

where, 

r - Optical path length in grams per square centimeter, 

j - Charge number, 0, . . . , J, where J is the charge number of the largest 
ion, 

Sj(E) - Stopping power associated with ion j, 

(j j - Total interaction cross section for ion j, 

0j(t, E) - Particle flux for ion j, 

Mj,* - Multiplicity, or the probability of creating fragments of ion j from an 
abrasive collision of ion A;, 

fj{ E) “ A known flux at the slab boundary as a function energy for ion j . 

The second term in the transport operator, ^S ; (E), is from the continuous 
slowing down approximation with the external force term in the original transport 
equation, equation (2.1). The external force being the electrostatic force between 
the target material’s protons and electrons and the ion’s protons. The stopping 
power of ion j is inside the differential, so a functional form of this term must be 
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found to continue the derivation. It can be related through a scaling constant to 
the stopping power of the proton, a well studied quantity, 

Sj(E) = FjSp(E), (2.3) 

where, V j is the scaling constant (Reference [19]), and Sp(E) is the stopping 
power of a proton as a function of energy. This approximation treats the incoming 
ion and target material as a sea of protons interacting separately instead of as 
discrete bundled atoms. The ion energies are large enough so that this 
approximation is valid. 

2.2.2. 1 Transformation from Energy to Path Length. When equation (2.3) is 
introduced into equation (2.2a), the stopping power of the proton is still inside the 
partial derivative term. This is transformed by scaling the energy variable in terms 
of the stopping power of the proton. The resultant variable transformation is from 
the energy of the particle to the energy integrated path length of the particle in 
the medium. The path length is defined as the penetration distance of a proton of 
initial energy Eo as it slows down to energy E 

S , S (E,E .) = (2.4) 

After changing the variable from energy to path length, the resultant set of 
equations is 

j 

Mj'k (T k rp k (T,s), 

*=7 + 1 


d _ d 
rr-‘'’Ts +a ’ 


(2.5a) 
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with boundary conditions of 

Vb(0,s) = (2.5b) 

The original energy dependent fluxes are obtained as follows 
V>j(r,s) = Sp(E) E), (2.5c) 

and 

f : (s) = Sp(E) fj(E). (2.5d) 

2. 2. 2. 2 Green's Function Solution to the Galactic Cosmic Ray 
Cascade Equations. General solutions to this set of equations for various 
boundary conditions can be found in Reference [6]; however, for clarity, a short 
synopsis of the Green’s function solution is given here. To generate an impulse 
forcing function, the boundary conditions are set to fj(s) = 6(s)Sjj . The Laplace 
transform of equations (2.5) is taken with respect to s and the resultant set of 
ordinary differential equations is solved analytically. The coefficients in the 
solution are dependent on the complex variable used in the transform. A recursion 
relation exists that can determine these coefficients (Reference [20]). Using a 
partial fraction expansion, the inverse Laplace transform can be determined 
analytically. In order to simplify the equation notation, all partial fraction 
coefficients are included even though some of the coefficients are identical or zero. 
Using an analytical mathematical argument, these coefficients can be identified. 
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The resultant equations determine the flux of any ion at any position and path 
length 

N, 


r=l 


(2.6) 


x nj ( s (« - "j- (, r ) - , 


where, 


N\ - The number of partial fraction terms for the J — l equation is 1(1 + 1)/ 2, 


" ^m(r) ^n(r)? 


* ^m(r) ^n(r)^ 

m(r) - Index: J — /, 

n(r) - Index: J — A; where = 0, 1, 
ii - Index: A: where k = 0, — 1, 

Z 2 - Index: /, 

- Modified partial fraction coefficients. 

For 1 < r < Ni-i, the recursion relation for the modified partial fraction 
coefficients are 


rJ—l 

*1 , r 


= Xu!,r E M j- l<k <r k Yf ur 9(Nj- k -r), 

jfc=j-i+i 


where, 


X 


j-t 


bj-U-u 

.Qj-U-^br — a T bj-u- 



(2.7a) 


(2.7b) 
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For Ni-i + 1 < r < fVj, the recursion relation is 

Hi = E (^r 1 ) Mi;J + c 

where, P r / selects only the non-zero values of With this Green’s function 

solution, all other distribution functions for fj(s) can be obtained with the 
appropriate integration over the Green’s function. 


2.2.3 Neutral Particle Transport 

The Boltzmann equation is simplified mathematically in Appendix 
Section A.l to describe the physical model for neutral particle transport. Steady 
state conditions, planer geometry, and isotropic scattering media are assumed. 

The source is an isotropically distributed source, and the particles cannot gain 
energy from scattering. The Boltzmann equation is then reduced to the 
heterogeneous slab, multigroup, neutral particle transport equation for group g , 
where g = 1, 2, . . . , G 

+ a ' 9 ^ X,fi ) = <t>g-(x,n') + ^s;(x), (2.9a) 

with a set of boundary conditions for each slab 

4> g {xi- ufi) = and 4> g {xi, — /x) = Fp l ’ , (/x), for /x > 0, (2.9b) 

where, 

i - Slab number, 1,2,..., NS, where NS is the total number of slabs, 


H - Angular direction = O • x, 
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^( x > A 4 ) - Steady state group flux in one dimension = /e * -1 dE $(r, fl, E, <), 

,• c , , , A /Jf s-, dE cr(r,E)$(r,fi,E,t) 

(7 - steady state group constant in one dimension = — — , 

/£-‘dE*(r,n,E,t) ’ 

<7*/_ a - Steady state group transfer cross section in one dimension, 


S^(x) - Steady state distributed group source in one dimension, 


Xi_i - Position of left boundary for slab z, 


x, - Position of right boundary for slab z, 


Fl'(aO ' Known general function of /z in positive half range (0 < /z < 1), 


FrXa 4 ) ‘ Known general function of fi in negative half range ( — 1 < fi < 0). 


2-2. 3.1 Singular Integral Equation Formulation. Equation (2.9a) can be 
solved in many different ways. The angular and spatial variables can be discretized 
and a finite difference formulation substituted for the derivative terms to create 
the Sn method. The angular dependance in the flux and source functions can be 
expanded in spherical harmonics to create the Pn method. Analytical methods can 
be used to solve for the flux. If the medium is infinite in extent in all dimensions, 
then a Fourier transform can be used to solve for the flux. If two semi-infinite 
mediums meet at the boundary, then a Laplace transform can be used to solve for 
the flux. The method of choice for this application is based on an integral form of 
equation (2.9a). To create a set of integral equations, equations (2.9) are 
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analytically continued into the complex plane. Through the Plemelj relations, the 
integrals axe then evaluated on the real axis. This is accomplished as shown below 
with details shown in Appendix Section A. 2 

<T* X 

Using equation (2.9a) \i is replaced with — /z and then multiplied by e 
The result is integrated over x on [zj, z 2 ] to obtain 


— -B g {n,s) - a' g dxe • <f> g (x,-fi) = 

3 J z\ 


= H - 7 T p ’ 3 A s ) + ^^-37 s p( 5 ^i^2), 


3 -=i 2 p ~ 3 


with, 


a g*2 


B g (fi,s) = e * 4>g(zi, -p) - e * <t> g (z 2 ,-p), 


PgAA = J Zi dxe f x d A <f>g'{x,-A), 

Ss(s,Z!,z 2 ) = f dxe ^S’(x). 


Next, fi is integrated on [— 1,+1] to obtain 

d P ~~ B A^A = AH 


+1 « p:as) 

r gg r V / i / \ , 

A g'g{s) + 


3'=1 


+ L(s)S'(s, Zj, z 2 ), 


where, 

AAA = Aa + 

A 

LM = i / +1 d/<— = ila(i±i), 

2 7-1 Li — S 2 \5 — 1 / 


Vff = 


i g' = g 
o a' # a 


(2.10a) 


(2.10b) 

(2.10c) 

(2.10d) 


(2.11a) 


(2.11b) 

(2.11c) 

(2. lid) 
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Equation (2.11a) is multiplied by e~ V 5- to obtain 


7h c ’ (l, ' s) = d £ ftK*. *»)*.-.(») + 

- 7 ' 1 ^ 3 g‘-\ 


where, 

Cg(n,s) = e • Bg{n,s) 


= - e - <t> g {z 2 ,-n). 


(2.12a) 


(2.12b) 


!«'(*> *i»*a) = 


°Vl 

e » 


Pa A*) 


1 /* + ! 

= - / dx e 
5 Jz, 


r+A 

J ipM 1 '—?)- 


(2.12c) 


To form a second equation, the signs of s and \i in equation (2.10a) are 


changed and the equation is multiplied by e to obtain 

J_ t „ _ -Dg(P’ s ) - a 'g Y! Jga'( S ’ 2 l’ Z 2)A i; / g (s) + 


3—1 
z 2 


+ L(s)S*(-s, 2 X ,z 2 )e * , 


where, 


(2.13a) 


D g (n,s) = -e B g (~n, -s) 

/ *q( z 2-* l) 


(2.13b) 


<T g 2 2 

Jgg f ( S '> Z U Z 2) = ™ / ? 53 '(' S ) 

1 *2 -*) f+l 

die 


1 r 2 3 *,(*»-*) r+i 

= - dxe z / d/i <£ fl ,(x,-/i). 

o */ Z] J — 1 

By defining j/q to be the positive root of the infinite medium dispersion 
relation, A gg (s), the integral equations, when evaluated at t/£, become a set of 


(2.13c) 



42 


constraints 


°) = (2.1 

L dfi J * h Q D ° M) = (2J 

The above integrals are not singular because if < 1 then |i/£| > 1 and real 
(Reference [21]). 

The integral terms in equations (2.12a) and (2.13a) are restricted to the 
real axis by use of the Plemelj relations (Reference [22]) 

.. 1 1 1 

um — ■ — — = p ± i7c6(ri — i/). (2.1J 

If the Plemelj relations are applied to this Cauchy type integral 


r+i 

X ( 5 ) = J drj 


Tj — S 


then, the result is 


lim *(i/ ± it) = X ± ( i/ ) = / + d t) ± /(„). (2. 

f — 0 J - 1 TJ — 1/ V 

All integration variables have been changed to 77 and i/ 6 [0, 1] U i / q . 

These rules are applied to the integral equations (2.12a) and (2.13a) to 


/" +1 n 

L d ” v) ± iirvC g (v,v) = 

= ^ E l 7A v ' z i’ z 2) + e*^S*(j/,2r 1 ,2 2 )L ± (77), 
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f + 1 77 

T <”/ - D g (r,,v) ± ixi/D g (i/, v) = 

«/ -~1 77 // 


~ 21 A +j(") + e ”^ _ S!;(-j/,Zi,Z2)L :t (j/), 


where, 


L ± (i/) = i / + ‘d, -i- ± " = ‘ 

2 7-i 77 — z/ 2 2 


1 , l—i/ ix 

2 n T+T * T ’ 




1 / ^ llTV(T 9'^9 

'W") ± - » , ■ 

2<7 <7 




i/<7*, r+l 1 

= V, + -£=* -f dr,-!— = 

lc rl 7-i 77 — 1/ 


c , , 1 - «/ 

S 9’9 + -itr In r— 
2(7* 1+1/ 


Eliminating I*j(i/, Zj, z-i) and J*^(i/, z lt z 2 ) by adding and subtracting the 
positive and negative branches of equations (2.16) to obtain 

j_ x ^ - zr - X aa{ v )Cg{v,u) = 

' a 9—9 

i ±i . 5 _, (2.1' 

5^5 a <? <;'=1 

/-I d '' ~ 5T S “ A »( i/ ) £, jK‘ / ) = 




3 S' = l 


Rewriting the integral equations in terms of the fluxes and changing the 
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integration variables so that they are evaluated on the interval [0, 1] gives 

i d, ^~M X " + 1 - 

A 

<t>g{Xi-l,-u) - <i> g (xi,-v)e~~ 


— e * 


- 

-1 L 


g^g g'=l 


i Zxfi-' 
a g e » 


V <7 


X{— i , X { ) , 


5— S 




7 

— e * 


1 rl 


2cr* 

“ T« ^99^) 


9^9 


<f>g{xi,v) - <f> g (xi„i,v)e~ 


r » 9 - 1 


S ^5'— ^i” 1 ’ * r *) 


n-t -t— ' 5 ’-*9 99 

9~*9 g f -\ 


*9 € * 


l/<7 


i ( £/, X,_i , X, ) . 


5—5 


(2.18a) 


(2.18b) 


where, Z\ and z 2 are defined at the slab boundaries, x,_i and x,-, and is defined 
as the dimensionless slab width a x g (x,* — x,_i ). A similar set of integral equations 
are generated for the interior slab points as shown in Appendix Section A. 4 


2*2.3. 2 The Fn Approximation. For the integral equations specified in the 
last section, various solution methods could be used to determine <^(x,/i). The 
solution method chosen for this treatment is a basis function expansion. At the 
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slab boundaries, the basis function expansions are defined as 

*»(*•-. -*0 = ^ eVv . M , 


and 


<7 cr=0 


r i Af-1 




5 o=0 


The basis function expansions for interior slab points are defined as 


and 


5 o=0 


r t jV-1 


where, 


J Or=0 


x, - Boundary slab points, 

Xj - Interior slab points, 
a&’* - Boundary expansion coefficients, 
b£'* - Boundary expansion coefficients, 
c 3 a ’ : - Interior expansion coefficients, 
d£ ,J - Interior expansion coefficients, 
V’o(i') - Basis functions to be selected. 


(2.19a) 


(2.19b) 


(2.19c) 


(2.19d) 
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When the basis function expansions are substituted into the singular 
integral equations, after much algebra, a set of equations is generated that specify 
the relations between the expansion coefficients and the various source terms for 


the slab boundary points 


N - 1 


a=0 L 


A* 


a *'B*W + b**A'(i/)e * 


— Bl* (//, x,_i ) + 


+ + — si;(z/,x,_ 1 ,x 1 ), 

0-0 & — 


1_ 

0-0 


(2.20a) 


iV-1 


Qr=0 


+ a^A5(i/)c- 


— R2* (z/, x,) -f 


1 


+ i T2 + • S2^(i/, , Xt). 


(2.20b) 
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For the interior, the integral equations become 


N-l 


C <y J ^ 9 a( L/ ) d£’ J A£(j/) — RV(l/,Xj) + — Tl^(l/,Xj,Xi) + 

a, — n /T * “ 


a=0 


N— 1 


3—3 


+ ^ b g -‘A g (*/), 

4^0 a=0 

E [ d a J B*(i/) - c"A*(*)] = R2'(i/,ij) + -E-T2‘(*.,x 1 -_ 1 ,x j ) + 


(2.20c) 
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a=0 


iV-1 


0—3 


+ -r-S2;(^,x,-„x J ) - <-.) 2 
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(2.21a) 

(2.21b) 
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(2.21c) 



47 


The inhomogeneous boundary terms that specify the source from adjacent 
slabs and the boundary conditions are 

= f dr, [Fft i ( 5 )C(»J(* j -x), ■,,») + 

i, . . 1 (2-22a 


^ *) = I'd* [Ff'fo) C(<r‘(l - t), „) + 

+ e-*-"' F’r (-7) S(<r;(i »)1 , 


where, 




S(-f , T7, i/) = 


e"«/" - 


7) — V 


TJ = 1/ 




1 — e 4 ■»+•' 

rj + v 


and x is defined as the slab boundary and interior points, or x £ {x,,x_,}. 

The down scatter terms that specify the source scattering from higher 
energy groups are 


X,X t ) — a g ^ Vg’—g x > x i ) i 


T2;(t/,x,-i,x) = (2.23b 

9 '= 1 

As shown in Appendix Section A. 3. 3, relationships for the integrals 
Ig'g'iViXjXi) and x,-i, x) are derived for two ranges of the angular variable. 
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When these integrals are singular, n = £ = s gg >£ < 1, then the relationships are 


„ i N-l 

~t -/ 


a 9* ->9^99'^' x ’ X *) “ 5 9 ^2 fc' l fca‘( S 99 , 0 "" X, X,) — 


cr=0 


s'-i 

” XT *V'-* x ’ Xt ’)’ 


^ N-l 


‘V— X i-X> X ) — 9 ; 9 ^ta( 5 W # 0 — ~S2^,(s^<f,X t _ 1 ,x) 


5 cr=0 


S"=l 


a a ® — ■£(’— 1 

^Or ^ ® ^ £« — If 


b® ’* x = £; 

d^ J l/l;. 


If s ss ,£ > 1, then the integrals are no longer singular, and the relationships 


(7 3 ^3 , 3'(' S 3S'£)I 33 '(£, x, £|) — (~ (s 33 '£)1 + 

a =0 

I pi* / „ C (^SS'OL^M'O 0, 1 / /• \ 

' K-VV^ss'S ) x ) — t ol 3 /(s 3fl -£, x, x, j — 

a 9 ’ 

g'-l 

®g ) ] ^fl' , 3 / (' S SS' < f X i 3<i)> 


a g^-9'9'( S 9a't)^gg>{ii X«-li x) — ( — s gg'i) + (■ s ps'<f )] + 

a =0 

+ R2j,(.„,f,*) - (3 ” ,01 ; (a ” ,0 S2^ (3 „.e.x,-„x)- 


° g XZ gg ft (^y X i-l •> X )t 

9 n - 1 


(2.25b) 
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where, 

1C 1 = 4 




— b^/e ~99 i 


X = X ,_1 


d« J - b a’‘e x) 


(2.25c) 


and 


^i = < 


J j , £ ( X * X < — 1 ) 

-a*’e 


X = X; 


a ; 7 o / » j { ( x ^'“l) / 

c a ~ a a e a * X ^ X,-. 


(2.25d) 


The terms that specify the distributed isotropic source for x,_i < x < x, are 


Sl^(i/,x,x,-) 


^ f X 'd z e~^ {z ~ x) S ‘ (z) ^ 0 

2/ «/ X 

Sj(x) i/ = 0, 


(2.26a) 


and 

( cr* /-a: <ri 

/ d2e-^ (r - z >SL(^) t/y^O 

17 - 71 ’- 1 (2.26b) 

s;(x) i/ = o. 

At this point, the source, S^(x), is a general function. Thus, the above integrals 
cannot be performed analytically in general. A numerical method of integration 
can be used to calculate values for this term. As an example, if a set of abscissas 
and source values is available that describes the source in the media, then an 
interpolation scheme coupled to a numerical integration routine can be constructed 
to evaluate the source terms. 


2. 2. 3. 3 The Post Processor Associated with the Fn Method. It has been 
shown by calculation that the Fn approximations, equations (2.19), converge 
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slowly with N . The convergence rate can be increased if the basis function 
summation term in the Fn approximations axe modified with the integral 
equations, equations (2.20). The first task is to regularize the singular B «(/i) term. 

The non-singular formulation of B® (7/) is determined from the original 
definition 

-7rr/d77— —Mv) 

g *'0 ^ 


■ (‘ ■ + ' £> ~) «'> - If £*, ^Mn). 


(2.27) 


If 


/ d t] — 

2 ( 7 ' J- 1 77 — 1/ 


'9 ' ,_1 *7 

is added to the second term and subtracted from the first term, then B£(t/) 
becomes 


B*» = 1+7/ 


g—3 

2cr‘ 


In 


1+7/ 


_ fj=. /‘ d , 

2(7' Jo T) — V 


(2.28) 


= A' (K)«.(r) - B?(K), 


where, A ffJ (i/) and are no longer singular at i/ = 77. L’Hospital’s rule is used 


to determine the integrand for B* 5 (z/) 




nj^aiv) - v^a{v) 

7] — V 


'M j T} = V 


(2.29) 


77 7^ 7/. 


The singular B ^(fi) term is eliminated in equations (2.20) by substituting in 
equation (2.28) to obtain 


7V-1 r 


E 

a=0 L 


{K.MM") - b;v)) + b;-'A»(^)e+ 


i._i) + — — + -i— Sl^K.ij.i.xiJ, 


9^9 


9^9 


(2.30a) 
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N-l 

E 

or=0 


K" {WM”) - B;*H) + a^Al(i/)e-^ 

= R2^(j/,Xi) + — —T2 t g (u,xi. 1 ,x i ) + E— S2^(i/,xj_i,x,). 


( 2 . 30 b) 


5^5 


9~*9 


These equations are solved for J2a=o a a' , ^’»(/ x ) an d I2a=o ^a'^aif 1 ) and 
substituted into the original Fn approximations, equations ( 2 . 19 ), to obtain 


</> 3 (xi_i,-/i) = F^'(/x)e + 

*b 2 <t' A* (/i) [ <T s- > s ^'^^ 7 x ‘ -1 ) "b 


( 2 . 31 a) 


+ Tijo., i,.„ *,) + £ {a}‘B;»(^) - b^AjMe-^ 


Or=0 k 


and 


M*,!*) = F£'*(/i)e-^ + 
1 


+ 


^p—oRS’ (/x, Xi) + S 2 ‘ (/r, x,-_t, Xj) 4 - 


2«i A «(M) 1 ^ ,vr ’ * 7 3 

+ T2;U,x,. 1 ,x,) + al_ g f: {b£* , 'B;'(/i)-a'- < A'( / i) 

o=0 l 


( 2 . 31 b) 


e “ 


The same procedure can be performed with the interior slab formulation to 


obtain 


= F6V)«‘ + , , . . 7 , x 

x x,) + + Tlj(fi, Xj,Xi) + 

+ d-»Ei<i J ' B 7 W + (<*;'- >) Ai( M ) 


( 2 . 31 c) 


a=0 


and 


+ rrrn t * 

2< n A ;» 

x [< 7 *^R 2 *(/i,x i ) + S 2 ;(/x,x i _ 1 ,x i ) + T 2 i(/i,x,-_ 1 ,x i ) + ( 2 - 31 d) 


iV-l 


a=0 


+ E dS J B;»(/i) + ctf - A'M 
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Equations (2.31) axe used to determine the final values of the angular flux. The 
derivation for the interior slab points can be found in Appendix Section A.5.1. 


2.2.3. 4 Extension of the Fn Method to Multiple Slabs. At this point, the 
incoming slab boundary flux values, Fr‘(/i) and F®’ *(/i), are general functions. To 
determine the fluxes for multiple slabs, the general functions Fr‘(/z) and F £’’(/•* ) 
must represent specific functions at the boundaries of the slabs. For the method 
being employed in this treatment, the rightmost incoming boundary flux, F^ NS (/r), 
is zero. The leftmost incoming boundary flux for group g is 

F£V) = S 3 o Hp - P a o), (2.32a) 


for a beam source, and 


Ft'M = So. 


(2.32b) 


for an isotropic source. The same source, either a beam or isotropic, is used for 
every group. To connect multiple slabs together, the i th slab’s incoming boundary 
flux is the outgoing boundary flux of the adjacent slab, or 

Fl‘(/*) = 1 ,/*) (2.33a) 


f r’(^) = -p) 


(2.33b) 


This connection is realized in the inhomogeneous boundary terms, 
equations (2.22). For a beam source incident on the leftmost face canted at the 
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angle / 4 > the terms axe 

RV g (n,x) = J^dr) rj [C(cr^(x,- — x), 77 , n)<f> g (xi, — 77 ) 


+ 


+ e * 


--*(x- 






+ 


(2.34a) 


, ,j C j -itL'/i+’il 1 - 1 '- 1 )) 5, i, ^ g , 

+ Ho S 0 e 0 S(cr g (xi - x),fi 9 0 ,n), 


and 


R 2 p(/^,X) = J 0 dr l *1 [ C (*g{x ~ X,-l),T/,/i)^(x,_i,J7) + 


+ e~^ x '- x) 


S ( &g ( 2- ^t-1 )l Tjl *?) 


+ 


(2.34b) 


+ a si - x^.^n). 

For an isotropic source incident on the leftmost face, the terms are 

Rl^,(/i,x) = j Q d V >? [C(<7^(x, - x),7/,/i)<£ 9 (x t ,-T?) + 

+ e -^(— *)s (<T «(x, -x^AUxi-i^) + 

+ So J^dT] A » +< '» (l_:r - l) ) S(< t‘(x,- - x),t?,/x), 


(2.35a) 


and 


R 2 ^(/i,x) = J dr] r/ [C(cr‘(x - x t _ 1 ), 77 ,/i)^(x i _ 1 , 77 ) + 

^(^g — g{ X ~ ^t-l)l * 7 ) l*)4 > g{ X i'> ~ T l) 


+ e 1 


+ 


(2.35b) 


a 

+ S o j Q dr l i7e“^*-» A » C(cr^(x - x,_ t ), 7, /x). 


The flux at the leftmost slab boundary, <^ 3 (xo, 77 ), represents the flux 


without the source, so 


<f>g( x o,v) = 0 . 


(2.36a) 
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Since there is no source on the rightmost slab boundary, the flux at that boundary 


is 


= 0 . 


These formulations are valid for slab boundary and interior points. 


(2.36b) 


2. 2. 3. 5 Scalar Flux Determination. Two methods are available for calculation 
of the scalar fluxes. The transport equation can be solved at fi — 0, or the angular 
flux can be integrated over \x. Both methods have approximations. The first 
method cannot use an exact value of zero for fx . A value of 1.0 x lO -20 is used as 
an approximation. The second method uses numerical integration for the basis 
functions. The direct integration method is the preferred method because it does 
not stress the post processor when calculating angular fluxes. The extra stress on 
the post processor, created when calculating the angular flux at a small //, could 
keep the Fn method from converging on an answer. 

The original transport equation can be used to determine the scalar fluxes 


^ 


<t>g{x,P) = 1 £ a 'g'^g / V <M x >/0 + xS‘(x). (2.37) 

a'= 1 J ~ l 

If fi is set to zero and the group scalar flux is defined as 

M x) = £V *,(*,#«'), (2.38) 

then the transport equations becomes 

<W*. o) = 5£4 - s <m*) + k(x). 

Z g' = 1 L 


(2.39) 



55 


Solving this equation for <f> g {x) obtains 

*.<*) = (2.40) 

9^9 fl g’=l a 9-+9 

This is implemented by first calculating the scalar flux for group one, then group 

two, etc However, the angular flux value for /x = 0, (f> g (x, 0), must be 

determined and that poses a problem as shown in equation (2.31). The 
denominator of some terms contain /x, so an approximation of setting /z close, but 
not equal to, zero is imposed and the resultant scalar flux value could be 
inaccurate or the Fn method could have problems converging while trying to 
calculate the angular flux at a small /x. 

The integration of the original Fn approximations, equations (2.19), over /x 
is the preferred method to obtain a scalar flux because most of the integrals are 
evaluated analytically. To obtain the scalar flux, Fl*(/x) and F j^ l (/x) are expressed 
in terms of known quantities at the slab boundaries 

fH/0 = Z] (2.41a) 

;=1 a g a =0 

and 

NS l N-l 

f r*(m) = Y “T E (2.41b) 

/=i+ 1 G 9 a=0 

These functions are substituted into the Fn approximations, and integrated over fi 
from [0, 1]. The result of this procedure is 

NS N-i 


r 1 ^ <r l N-l 

4>g{x 0 ) = j dr/ U(tj) + Y^f 1 Y &3 a x 

J0 (=1 a 3 cr=0 

x J dr/ 0 a (r/)e - ^ A », 


(2.42a) 
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/•l NS t NS -l N— 1 

<f>g(^Ns) = I dr} U(r})e~* ^ ^ bj* x 

J ° 1=1 ^3 o =0 

x J^dr] xp a (T])e~n ■ 


(2.42b) 


_iv NS A * 


M*i) = A lW(,)e-'C..^ + b o*' + ^a^ +1 l x 

70 o^o L °*g a 7 \ 

x / dr/i/> a (r/) + J2 b 9 a '‘ f dr/ 0 Q (r/)e“* £*='+i A £ + 

Jo 1=1 a g o=0 Jo 

NS l N- 1 x 

+ II - J= j 3 ‘ E / d7/t/> a (77)e"^*=-+i A 5, 

/=»+ 2 ^3 o =0 - 70 


and 


^s( x i) = / dr/ ZV(r/)e aE*=i A * + *(->)] + 

er* *-i . .-1 / AT-i 

-2=2 2 (d« + c") / d, «,) + £ -*7* £ X 

a o=o - 70 ;=i a g o=o 

x J dr/ r/r a (r7)e - ^E*=‘+ I + <T s>( r .' -Xi - 1 )] -j. 

NS N- l I _ l 

+ E-?Ei: J / d, <M,)e~ P-- 

(=,+l <*3 o =0 ' / ° 


A 5 + 


where, 


U{tj) = ^ 


(2.42c) 


(2.42d) 


So S(t) — Hq) Beam Source 
So Isotropic Source. 

All integrations over the basis functions are performed numerically using the 
standard shifted Gauss-Legendre quadrature scheme. 


(2.42e) 


2.2.4 Coupling of the F N and Galactic Cosmic Ray Cascade Models 

From the physical model, a distributed monodirectional source due to the 
galactic cosmic ray cascade exists within the slab. The F N algorithm, as derived, 
allows for a distributed isotropic source. Since a direct substitution of the cascade 
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source into the Fn algorithm cannot be performed, generation of an isotropic 
source from the monodirectional cascade source is used as an alternative. This is 
accomplished by determining the collided flux from the cascade source and setting 
the Fn source to be this collided flux. This derivation is performed for a single slab 
only. 


2.2.4. 1 Generation of the Fn Isotropic Source. The formulation of the 
isotropic Fn source term in equations (2.9) starts with the uncollided and collided 

solutions to the transport equation with a monodirectional source (0 < //git < 1) 

d I s y+i 

Q + <Jg (j)g(x,n) = — Gg'—g I d H <f)gl(x, (i ) + 

J *'-1 (2.43a) 


+ ~ ^git), 

with boundary conditions of 


<t> g (x o,//) = 0 (2.43b) 

where, ^ s (x) is the neutron flux found from the galactic cosmic ray cascade model 
as described in Section 2. 2. 4. 2. 

Let the angular flux be composed of the uncollided and collided angular 

fluxes 


<j>g(x,n) = <t> u g {x,n) + <t> c (x, n). 


Substituting into equations (2.43) yields 
d 

* Tx + 


(2.44) 


<t>g(x,P) = Vg(x)6(fi- GIT), 


(2.45a) 
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<f>g{XO,H) = 0 , 


(2.45b) 


and 


»Tx + '' 


<f> C g(x,l*) = I Z <V-» / + 

S'=l - 7-1 

+ ^ Z <v-» / V <P}(x,n'). 

S'=l J ~ 1 


(2.45c) 


Equation (2.45c) is the transport equation used in the Fn method for a single slab 
if 


S s( x ) = Z a 9'-9 I <V 4> u g ,(x,n'). 

g‘=\ J ~ X 


(2.46) 


Therefore, to obtain the Fn source term, equation (2.45a) is solved for the 
uncollided flux, which is then substituted into equation (2.46). 

A solution to equation (2.45a) for a single slab is 


^ g 

= f dz- 'i> g (z) 6 (fi - hgit)- 

Jxo /•* 


Substituting equation (2.47) into equation (2.46) yields 

S s( x ) = Z 'V-s I <V [ dz - ^ g ,(z)8(n' - ficn). 

g'-l J ~ 1 Jl 0 /* 

The integral on fi is easily performed to obtain 


(2.47) 


(2.48) 


S f (*) = 


— Z a 9'~9 / 

g‘—\ Jx 


v ( *~° 

: e **GIT 

dz * a .(z). 

*o A*GIT 


(2.49) 


This formulation can then be placed into equations (2.26) to determine terms 
required for the Fn method for a single slab. The next task is to determine ^(x). 
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2.2. 4. 2 The Galactic Cosmic Ray Cascade Formulation for the Fn Source. 

The flux found from the galactic cosmic ray cascade transport equation, 
equation (2.5), is the flux of particle j based on the particle fragmentation and 
transport through the medium. The flux that must be used as the Fn source is for 
neutrons which is found by setting j to 0. Because the neutrons are acting as a 
source, the transport part of equation (2.5) is set to zero so that particles are not 
transported in space and energy by the galactic cosmic ray equations when 
created. The resultant algebraic transport equation is solved for the neutron flux 

1 J 

Vg( T ) = ipo(T,E g ) = — <7k EJ, (2.50) 

*0 k=2 


where, E g represents the energy of the group being calculated. 

When the general Green’s function solution, equation (2.6), is introduced 
into the above equation, the result is 

^s( T ) = 


i J—1 _ N x 

i v . - CTj-l 


bp(t< s ) (=1 < 7 q r=1 


(2.51) 


* Yf,"' [d(j, - - e(s, - l , 

where, s g is the path length evaluated at E g . 

To convert from the optical path length, r, to a physical space dimension, 
z , this differential relation must be satisfied 


*.(T)|dr| = ¥ s (z)|d* 


If r = pz, where p is the density of the target material, then 


(2.52) 


V g (z) = pV g (T). 


(2.53) 
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To complete the variable change, these definitions must be specified 


(Tj-l = pcrj-i. 


Sp(E 3 ) = /?Sp(E s ), 


y*E° f 1 l 

~ h a dE pSp(E 3 -) = ~p Sa ' 


Therefore, the final expression for the galactic cosmic ray cascade source used in 
the Fn source is 

(=1 r=1 , 9 ... 
r 1 l“’ 0 ' 

^ \rJ-l I at _ S 9 _\ a< S 3 _\ 


This result is used in equation (2.49) to generate a source of neutrons from the 
galactic cosmic ray cascade. 


2. 2. 4. 3 Fn Distributed Source. In the last two sections, an Fn neutron source 
is derived from the galactic cosmic ray cascade. This section places the result in 
the Fn context. 

The cascade neutron source in equation (2.55) is substituted into the Fn 
source term in equation (2.49) to obtain 


s,(*) = — EMo ^-lEYfc'x 

^0/iGIT , =1 r=1 


X! c /I? \ *1 ? *0 2 > ^)] > 

g '= 1 


X 
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where, 


I 9 '(Z,x) = /d ze _*). 

Jx 0 € 


(2.56b) 


To complete the analysis, equations (2.56) are substituted into the Fn 


source term, equations (2.26), for a single slab to obtain 

J—2 Ni 


SI (j/,x,xj) = 


P°g 


VqPGYIV (_! 


J 2 £ X 


yj-l y' g ~*g 

* i,r h me,) 


ipi, 


S2 (z/,x 0 ,x) = 


P°9 


J—2 


r = l 


Ni 


(2.57a) 


C’o/iGIT*' ;=1 


m 0 t j-i<Tj-i y x 


r = l 


X 


<7g'—g 


(2.57b) 


y j - 1 y 

n ' r £Sp(E s -) 


IP2, 


where, 


IPI = r'dze-T?' ! -l[V(Pj-n,2)-V(^ 1 .*)l, 

J X 


(2.57c) 


and 


IP2 = rdze-^*-*[I gl {Vj- il ,z)-I g .{Vj^z)\. (2.57d) 

J Xo 

For v = 0, the term in equations (2.26) is the source value in equations (2.56). 

The integrals in equations (2.56b), and (2.57) have been evaluated using 
the Maple symbolic manipulation routines in MaTHCAD (Reference [23]). 
However, the integrals in equations (2.57) must be partially evaluated to allow 
evaluation in MaTHCAD. The full integrals are 


IPI 


= ndze -^’-* ) /d 

J X J Xn 


. v J-< 


w -W^- 


. V J-<1 


— W 


(2.58a) 
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and 


IP2 


= [*dz [ Z (lw e ^git (z u ' ) ~ vV x 

</rn An 


* Xq J Xo 

x [* f - uA - 0 ( V 


— w 


L V^-i ) \ v J-i2 J 

The partially evaluated integrals are 


IP1 = - r i d2e‘^ (j ‘ l) '^ ,z ”' ) 

J q\ J x 


- [ n dw [ X 'dz 

J q 2 J W 


Xl j — r) 2 — (z — o') 

y| — V ' 1 Lkr* T'-n ' ' 


IT 


and 


IP2 = 


-r 

Jq x 


du; e-^-M )"- 


■f?V f 

J W 


j — ^(r— z) 2 — (z— it/) 

dz e u ^git , 


where, 


px = min(^^ — ,x) 
^J-» 2 


<7x = max(^ — ,x 0 ), 
v J-h 


and 


p 2 = min( _ 9 , x x ) 

Vj-i^ 


V 


q 2 = max( — , x). 

"J-ii 


These are the integrals symbolically evaluated by MaTHCAD. 


(2.58b) 


(2.59a) 


(2.59b) 


(2.59c) 


(2.59d) 
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Chapter 3 

Numerical Implementation 

The numerical methods required to solve the Fn equations are simple and 
straight-forward. The tasks that require some form of numerical method are 
determination of the expansion coefficients, collocation points, Legendre 
polynomials, general source function, basis function, and various integrals. The 
following numerical algorithms are required: matrix inversion, cubic spline 
interpolation, and Gauss-Legendre integration. Also, two iterations are required to 
complete the algorithm. 

3.1 Evaluation of Expansion Coefficients 
The integral equations, equations (2.18), are in the form of the 
inhomogeneous Fredholm equation 

Hx) — A / d£ k(x, £)<!>(£) - f(x) = r(x). (3.1) 

J a 

If the solution, 4>(x), is replaced with a basis function expansion 

n 

<f>(x) = ^ajUjix), (3.2) 

j=o 

then the resultant equation is 

- x it a i / d £ - f(x) = r (x), 

J=0 j = 0 Ja 


(3.3) 
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where, the kernel, fc(x,f), is singular but can be regularized. 

To determine the coefficients, equation (3.3) is multiplied by a set of test 
functions u;(x), i = 0, 1, . . . , N. The inner product of v,-(x) and r(x) is then 
required to be zero for each value of i 

f d xv{(x)r(x) = 0 i = 0, 1, . . . , n. (3.4) 

Jo 

If U{(x) is equal to 5(x — x,), then the process is called the collocation method . If 
Vi(x) is equal to x\ then the process is called the method of moments. A more 
general method called the Galerkin method uses a general weight for t\(x). The 
collocation method was chosen for this problem because the above integrals can be 
performed analytically. The matrix elements are the integrals with the rows 
representing the order of the basis function expansion j and the columns 
representing the basis of the delta functions i. The solution vector is the basis 
function expansion coefficients. The inhomogeneous vector is composed of the 
terms associated with the boundary conditions, external sources, down scatter 
sources, and slab interfaces. 

From experience, the choice of the collocation points, x, in the above 
formulation and vq in the F n notation, has not greatly affected the solution or its 
rate of convergence, as long as the positive root of the infinite medium dispersion 
relation or constraint, i/q, is included. Three sets of collocation points have been 
used extensively in Fn calculations. The first set is N — l evenly spaced points 
over the interval [0, 1]. The second set is formed by the zeros of the Legendre 
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polynomials up G {Pn~i(2u — 1) = 0}. These are determined from a subprogram 
in Reference [24]. The third set is formed by the zeros of the Chebyshev 
polynomials up € {T^-\{2u — 1) = 0}. These are determined from 


1 1 

v <‘ = 2 + 2 C ° S 


2/3 — 1 
21 V : 


(3.5) 


Recall that i/q is the positive real root of the infinite medium dispersion 


relation 


A flP K ) = 0 = 1-4 ^ f ln 


i-4 


i + 4 


(3.6) 


The behavior of this function is well known. In the interval [l,oo), the function 
monotonically decreases; therefore, there is one real root. The function is 
symmetric about zero, so it has another real root in the interval (— 00 , —1]. To 
determine the positive root, a routine is used that is based on the bisection root 
finding technique. The condition number of the collocation matrices is very 
sensitive to the precession of v q, thus, a relative error equal to the machine 
precession is required. This appears to give acceptable results over all ranges of 
operation. 

The use of the collocation method transforms equations (2.20) into two 


matrix equations of the form 
N - 1 

/* ^ \p‘cx'—>a(3 ^a^Qr/ 3 ] = T^, 

Q=0 


(3.7a) 


and 


N - 1 

^ ^ [^or^or/3 “I" a@\ — r^, 

Or=0 


(3.7b) 
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where, 

Vc .0 = A*(ty), (3.7c) 

and 

=- a0 = BIM. (3.7d) 

These matrix equations could be solved directly, but tya/? goes to zero as a 
increases. This creates an ill-conditioned matrix equation. A change in the form of 
the matrix equations not only side steps the ill-conditioned matrix problem, but 
reduces the computer storage requirements by a factor of two. To reformulate the 
matrix equations, equation (3.7a) and equation (3.7b) are added to obtain 

N - 1 

52 [(o a + b a ) (~ a 0 + ^ 0 /j)] = (T^ + r^). (3.8a) 

o=0 

The same equations are subtracted to obtain 

N-l 

52 [( a Q ~ b a ) (z. a p - 'Jo./?)] = (T/j — T^). (3.8b) 

o=0 

The new matrix equations are solved for g£ and g~ by letting 

9a = + (3.9a) 

and 

g~ = a a - b a . (3.9b) 

The original coefficients can be recovered from 



(3.10a) 
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and 

b ° = \ (at - 9«) • (3.10b) 

The matrix equations are solved using LU decomposition as found in 
Reference [24], which is based on C rout’s algorithm with partial pivoting. The 
coefficients are determined using the companion back-substitution algorithm. The 
matrices only need to be decomposed when the number of terms used in the basis 
function expansion changes. 

The QR and Singular Value Decomposition (SVD) algorithms from 
Reference [24] were used instead of the LU decomposition algorithm when Fn 
algorithm did not converge; however, the QR and SVD algorithms did not make 
the Fn algorithm converge either. Therefore, the matrices created in this 
algorithm are not ill-conditioned. Since the LU decomposition algorithm is a fast 
and convenient algorithm, it is the algorithm of choice. 

3.2 General Function Integration 

The terms containing integrals that must be evaluated either analytically or 
numerically are the matrix terms A£(z/), B£(z/), and B* 5 (^), the inhomogeneous 
terms Rl^(j/, x) and R2 l g (v, x), and the source terms Sl^(i/, x, x t ) and 
S2*(i/, x,_i , x). Due to the generality of the functions involved, the Gauss-Legendre 
quadrature is used. The quadrature points and weights are found with the GAULEG 
subprogram from Reference [24]. 
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An alternative method to determine A^(j/) and B£(j/) can be found using 
the Legendre polynomial as the basis function. The Legendre polynomial recursion 
relation is used to determine a recursion relation for these terms. Therefore, 
numerical integrations need not be performed. The derivation for these relations 
are found in Reference [12]. 

The recursion relation for A£(£) with a > 1 is 

a;({) = )a;_ ! (o + 

j V a J \ a J (3.11a) 

and 


AS(f) = 1-fln 



(3.11b) 


The forward recursion is stable only for £ 6 [—1,0], and the backward recursion is 
stable for £ ^ [— 1,0] which can be evaluated with Miller’s algorithm 
(Reference [24]). 

The recursion relation for B® (f) with a > 1 is 


£>;«) = (^-i) t 2 ? - !) - (^-) b ;. 2 ( o+ 


9~+9 £ ” 9^9 T 

°a,2 

rr* 


2(j. 


and 


(3.12a) 


BS(f) = 2- 


9-^9 


1 +£ln 

1 + 7 

. 

{ J 


(3.12b) 


The forward recursion is stable only for £ £ [0, 1], and the backward recursion with 
Miller’s algorithm must be used when £ = i/q. 
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To evaluate B* s (^), equation (2.28) gives 

- B»», (3.13) 

where, B£(/i) and 'ifc ot (n) are determined from the appropriate recurrence relations. 

3.3 Determination of the Legendre Polynomials 
The standard recursion relation is used to determine the value of the 
Legendre polynomials for any value of x. This function is used for determination 
of the basis functions 

Pn(x ) = x— -P^x) - — P n _ 2 (x), (3.14a) 

Tt Tt 

Po{x) — 1 and Pi(x) = x. (3.14b) 

This formulation provides a stable calculation of the Legendre polynomials over all 
orders, n, and all values of x being used in the Fn algorithm. 

3.4 Cubic Spline Function Interpolation 
The Fn formulation contains a general isotropic source distribution as 
shown in equations (2.26). This source is specified with two vectors for each 
energy group with a source. The first vector contains the slab positions of the 
source values. The second vector contains the source value at the slab positions. 
Since this source must be integrated at various Gauss- Legendre quadrature points, 
a natural cubic spline is fit to the source values. The algorithms that determine 
and store the second derivative values and interpolate the source values are from 
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Reference [24]. The Gauss- Legendre integration algorithm for equations (2.26) 
uses pre-calculated source values interpolated at the quadrature points. 

3.5 Basis Functions 

In contrast to the choice of collocation points, the choice of the basis 
functions, 0 O (^), is crucial to the convergence of the Fn method. The use of 
monomials, /*“ and (2{i y — l) a , where 7 is a value between 0 and 1 , has been 
shown to generate ill-conditioned matrices (condition number ||A||oo >> 1), and 
therefore, unstable convergence properties. The use of shifted Legendre 
polynomials, P a ( 2 /i — 1 ), has been shown to have fairly stable convergence 
properties. However, modified Legendre polynomials, P a ( 2/z 7 — 1 ), when they 
converge, tend to converge faster than the shifted Legendre polynomials. 

The derivative of the basis functions is required in equation (2.29). The 
table below shows some examples of basis functions and their corresponding 
derivatives. 


3.6 Iteration Schemes 

Two iterations are required to facilitate the Fn algorithm. The first 
(Iteration I) is an iteration to converge the fluxes over the number of basis 
expansion terms, N. The second (Iteration 2 ) adjusts the interior slab boundary 
fluxes with constant N to accelerate the first iteration. 
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3.6.1 Iteration 1 

As the number of expansion terms is increased, the relative error between 
the fluxes calculated from the current and previous expansion terms is determined 
for each boundary point and angle. If the maximum relative error is smaller than a 
prescribed value, then the boundary problem is deemed converged. If the 
maximum relative error is larger than the prescribed value, the number of 
expansion terms is increased, and new fluxes are calculated. 

Since the interior integral equations depend on the boundary expansion 
coefficients, the following algorithm is imposed. All slab boundary expansion 
coefficients must be determined for every iteration in N. If the boundary fluxes do 
not converged, then computer time is not wasted computing interior fluxes with 
inaccurate boundary coefficients. Once the boundary fluxes converge, then all the 
fluxes are determined, and the convergence of all points is tested as described 
above. 

3.6.2 Iteration 2 

To solve for the fluxes in slab i, the right boundary incoming flux from slab 
i + 1, <f> g (xiy — /i), is needed. But it has not yet been determined. Therefore the i tK 
slab’s boundary fluxes are inaccurate. If multiple calculational passes are made 
through the slabs, then the inaccuracies are diminished. There are three methods 
that could be used to perform these passes. The first is to create a slab based 
matrix and use a solution matrix iteration technique, such as Gauss-Sidel, or a 
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simultaneous technique, such as LU decomposition, to solve for the resultant 
boundary fluxes. The next method determines the fluxes using multiple passes 
until the difference between successive passes is less than a prescribed tolerance 
value. The final method calculates the fluxes for a specified number of passes. 

The first and second methods obtain the best answer for the chosen number 
of expansion terms. However, in the iteration over expansion terms, the most 
accurate interior boundary fluxes are not necessary. A balance must be obtained 
between the time of calculation of the most accurate slab boundary fluxes for a 
certain number of expansion terms and the determination of the number of 
expansion terms that converge the boundary fluxes. The third method achieves a 
sufficient balance by estimating the interior slab boundary fluxes and allows the 
iteration over the number of expansion terms to drive the final boundary fluxes. 

With these numerical methods and considerations, a computer program has 
been written to implement the Fn algorithm. In Chapter 4, the program is verified 
by performing internal consistency checks and comparisons to standard codes. 



Table 3.1: List of Basis Functions and Derivatives 


Number 

'l’a(v) 


i 

P.{2u- 1) 

!_ ( £L 1)a - (21/ - l)0 a (l/)] 

2 

Pa( W - 1) 

[fe-.M - (2W - 1 )».(>)] 

3 


O0o-1 (v) 

4 

(2i/' r - l) a 

2i/ 7-1 a7^ Q -i(i/) 
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Chapter 4 

Program Verification 

Three separate programs are available for application of the Fn algorithm. 
The MGSLAB program is for heterogeneous media and a single slab for the galactic 
cosmic ray cascade source. The MGSEMI program is for a homogeneous 
semi-infinite media. The FNCRIT program, for verification only, determines the 
critical slab thickness and fluxes for a single slab and energy group. Even though 
the FNCRIT program is not applied to the galactic cosmic ray cascade problem, it 
is included in the Fn algorithm package created. These programs must be verified 
internally and against other verified and accepted programs. Section 4.1 describes 
the internal verifications performed and associated results. Section 4.2 describes 
the external verifications performed and associated results. 

In this chapter, it will be shown that the MGSLAB, MGSEMI, and FNCRIT 
programs are internally consistent and any discrepancies between these programs 
and ANISN/Pc (the program used for external verification) are well understood 
and fully explained. Therefore, the programs are viable tools that can be used to 
analyze the galactic cosmic ray cascade and other problems. 
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4-1 Internal Verification 

Various measures are used to verify that the heterogeneous finite and 
homogeneous semi-infinite mathematical models are programmed properly in 
MGSLAB and Mgsemi. FNCRIT is not internally verified because it is the 
MGSLAB program with modifications necessary for inclusion of the critical slab 
calculation. The first and most important internal verification is the proper 
compilation of the programs. In this and all other sections, it is assumed that the 
programs are compiled and linked properly. 

Six measures are used to verify proper programming. The first measure is 
verification of faster convergence rates and more accurate fluxes for increasingly 
accurate computed integral terms. When the integrals in equations (2.20) are 
computed with higher accuracy, the Fn algorithm should converge with fewer 
expansion terms and the fluxes calculated should be more precise. The next 
measure is more accurate and consistent fluxes for a decrease in convergence 
tolerance. As the precision of the converged fluxes is increased, by decreasing the 
convergence tolerance, the fluxes should approach the same value. The third 
measure is a consistent number of expansion terms for convergence with variation 
in the initial number of expansion terms. The problem should converge with the 
same number of expansion terms without regard to the initial number of expansion 
terms. The next measure is consistent fluxes with variation in basis functions. The 
resultant fluxes should be the same no matter what basis function are used. The 
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next measure is verification of a faster convergence rate for an increase in the 
number of passes for the internal slab boundaries. As the internal slab boundary 
fluxes approach their true values, the number of expansion terms for convergence 
should decrease. The last measure is consistent fluxes when slab boundary flux 
calculations are compared to slab interior flux edits. The algorithm should 
calculate the same fluxes whether the interior or boundary formulation is used. 

4.1.1 Variation of the Number of Points Used in Quadrature Integrations 

As the number of quadrature points is increased for the matrix element 
integral evaluations, the fluxes should become more accurate, and the Fn algorithm 
should converge faster because fewer expansion terms are required. This is 
investigated by increasing the number of quadratures points for the matrix element 
integrals, and noting the number of expansion terms required for convergence. 

The physical problem being analyzed for MGSLAB and MGSEMI is a single 
energy group problem with a beam source of unit strength normally incident on 
the left face of the slab. For MGSLAB, a single slab one centimeter thick is used. 
The total cross section is 1.0 cm -1 and the scattering cross section is 0.99 
cm — steradian -1 . Shifted Legendre polynomials, P a (2/i — 1), are the basis 
functions. The convergence tolerance is 10~ 6 . 

The number of quadrature points is increased from five to ninety as shown 
in Tables 4.1 and 4.2 for the MGSLAB and MGSEMI programs, respectively. From 
these tables, as the integral evaluation becomes more accurate, fewer expansion 
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terms axe required for convergence of the problem and the fluxes become more 
accurate. 

4.1.2 Variation of the Convergence Tolerance 

The iteration scheme which determines the final number of expansion terms 
used in the Fn algorithm is controlled by the convergence tolerance. A relative 
flux is calculated between successive iterations and compared to this tolerance at 
every edited spatial and angular point. As this tolerance is decreased, the fluxes 
should approach the same value. This is investigated by decreasing the 
convergence tolerance and comparing fluxes at three different spatial points and 
one angular point. 

The same physical model is analyzed as above. The number of quadrature 
points used for the matrix elements integrals is thirty-five. The reported fluxes are 
converted to units of flux per steradian. 

The convergence tolerance is decreased from 10~ 2 to 2 x 10“ 8 and compared 
at fj. = 0.5 for x = 0.25, 0.50, and 0.75 centimeters as shown in Tables 4.3 and 4.4 
for the MGSLAB and MGSEMI programs, respectively. From these tables, as the 
convergence tolerance is decreased, the fluxes converge to the same number. 

4.1.3 Variation of the Initial Number of Expansion Terms 

The iteration scheme which determines the final number of expansion terms 
used in the Fn algorithm must start with an initial number of expansion terms. 
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For a single slab, the convergence properties should not be a function of the 
number of terms used to start the iteration. This is investigated by varying the 
initial number of expansion terms and noting the number of terms required for 
convergence. 

Using the same physical model as described in the sections above, the 
initial number of expansion terms is increased from three to twenty-one as shown 
in Tables 4.5 and 4.6 for the MGSLAB and MGSEMI programs, respectively. From 
these tables, as the initial number of expansion terms is varied, the number of 
terms required for problem convergence does not vary. 

4.1.4 Variation of the Basis Functions 

As the basis functions are changed, the fluxes should remain the same. This 
is investigated by varying the basis functions and comparing the resultant fluxes. 
The basis functions being varied are the modified Legendre polynomials, 

P a (2fj,' y — 1), where 7 is a value between 0 and 1. The same physical problem 
above is being analyzed. 

The adjustable parameter, 7, is varied from 0.70 to 1.0 for MGSLAB and 
0.65 to 1.0 for MGSEMI as shown in Tables 4.7 and 4.8 for the MGSLAB and 
MGSEMI programs, respectively. From these tables, the fluxes are the same no 
matter what basis function is used. The number of expansion terms required to 
converge the problem however, varied depending on the basis function used. This 


is an expected result. 
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4.1.5 Variation of the Number of Inner Slab Boundary Iterations 

For multiple slabs, an inner iteration which estimates the interior slab 
interface fluxes is used to accelerate convergence in N. With each execution of the 
Fn algorithm over the slabs without changing the number of expansion terms, the 
interior slab interface fluxes approach the true fluxes. Therefore, as the number of 
inner iterations is increased, the Fn algorithm should converge faster by requiring 
fewer number of expansion terms for problem convergence. This is investigated by 
increasing the number of inner iterations and noting the number of expansion 
terms required for convergence. 

The physical problem being analyzed for MGSLAB is a single energy group 
problem with a beam source of unit strength normally incident on the left media 
face. Four, one centimeter thick slabs of the same material used in the previous 
verifications is specified. The convergence tolerance is 10" 4 . 

The number inner iterations is increased from one to five as shown in 
Table 4.9. From this table, as the number of inner iterations is increased, fewer 
expansion terms are required for problem convergence. 

4.1.6 Variation of the Number of Slabs and Interior Points 

Whether the Fn algorithm determines the fluxes at interior slab (edit) 
points or at slab interfaces, the final flux value should be the same. This is 
investigated by executing MGSLAB with one energy group and four slabs of the 
same material with no interior points and with one slab with three interior points. 
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The same physical problem analyzed above is used. 

The results are shown in Table 4.10. From this table, the flux values axe 
identical whether calculated as slab interface points or interior slab points. 

4.2 External Verifications 

In the last section, it was shown that the mathematical model is 
programmed properly. This section shows some of the verifications used to 
compare the Fn algorithm to Anisn/Pc (Reference [25]) which is a one 
dimensional, anisotropic, multigroup, Sn transport code for the Pc computing 
environment. ANISN/Pc has been modified to execute on a Vax/Vms machine, a 
Dec ALPHA machine with the OPEN VMS operating system, and all UNIX 
machines. 

4.2.1 Finite Slab Comparison 

The MGSLAB and ANISN/Pc programs are used to determine the fluxes for 
all combinations of one to three energy groups and one to three slabs for an 
incident beam of unit strength and an exponential distributed source of the form 
e~ r . The exponentially source is fit with a cubic spline and the source term 
integrals evaluated with a Gauss-Legendre quadrature order of fifty-six. For all 
cases, the transport media is one centimeter thick. For two slabs, the slab 
thicknesses are 0.4 cm and 0.6 cm. For three slabs, the slab thicknesses are 0.2 cm, 
0.4 cm, and 0.4 cm. The convergence tolerance is 10~ 4 . The number of quadrature 
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points for the matrix elements and the boundary terms is thirty-five. The 
Anisn/Pc program is executed with sixteen discrete angles and 100 mesh cells. 
The number of direction edit points in MGSLAB are thirty-two and the number of 
slab position edit points is forty for a single slab, fifty for two slabs, and sixty for 
three slabs. The beam source incident on the left side of the transport media is 
directed along the closest discrete angle represented in Anisn/Pc to the normal at 
the slab face. A slim, 10 -6 cm, vacuum slab is constructed in the ANISN /Pc 
model to contain the shell source which models the delta function source used in 
Mgslab. 

The cross section sets used are consistent for all cases. For the first region 
or slab, the cross sections are 


g 

•i 



a 3~r 

1 

2.0 

1.0 

0.0 

0.0 

2 

5.0 

1.0 

3.0 

0.0 

3 

3.0 

0.1 

0.2 

1.5 


For the second region or slab, the cross sections axe 
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To transform the angular fluxes into consistent units of particles per square 
centimeter per second per steradian for comparison, the ANISN/Pc values are 
divided by 4 tt, and the Fn values are divided by 2 tt. 

The comparison between Anisn/Pc and MGSLAB angular and scalar fluxes 
for three groups and three slabs with a beam source are shown as an example in 
Figures 4.1 and 4.2. Figures 4.3 and 4.4 show the example comparison plots for 
the distributed source. The lines represent the MGSLAB fluxes, and the symbols 
are the ANISN/Pc fluxes. For each group, the slab boundary fluxes at 0.0 cm, 

0.2 cm, 0.6 cm, and 1.0 cm along with the fluxes at the slab centers of 0.1 cm, 

0.4 cm, and 0.8 cm are shown on the angular flux plots. The MGSLAB angular 
fluxes match the Anisn/Pc angular fluxes except along the direction of the beam 
source. This is because Anisn/Pc prints the collided and uncollided angular 
fluxes added together, whereas, theoretically, the uncollided flux along the beam at 
the left slab face is infinite and exponentially decreases inside the slab. This effects 
the scalar flux in group one because ANISN/Pc is trying to model the delta 
function source. Figure 4.2 shows this discrepancy. When a distributed source is 
used, as in Figure 4.4, this problem is not encountered. 

From the example plots shown and those not shown, MGSLAB and 
Anisn/Pc generate the same fluxes for the same physical problem. This verifies 
that MGSLAB can be used as a tool for solving heterogeneous slab neutral particle 


transport problems. 
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4.2.2 Semi-infinite Media Comparisons 

The MGSEMI and ANISN/Pc programs are used to determine the fluxes for 
a semi-infinite slab with the same beam source used in the finite slab comparison. 
Only one region is used, but one to three energy groups are compared. The same 
cross sections are used as in the finite slab case. The ANISN/Pc and MGSEMI 
angular flux values are compared only at the slab boundary. The scalar flux values 
are not compared. The limited comparison is due to the ANISN/Pc program 
generating results that are influenced by the rightmost boundary and the beam 
source. To overcome this limitation, the slab thickness in ANISN/Pc is greatly 
increased. Figure 4.5 shows an example comparison plot for three groups. The 
group one fluxes are not identical, but have the same shape. The group two and 
three fluxes are approximately the same. The discrepancy in group one is mainly 
due to the beam source delta function problems as seen in the finite slab case. The 
delta function is smeared out over all angles because of the infinite extent of the 
slab for the scattering process, unlike the finite slab case. 

From the example plots shown and those not shown, MGSEMI and 
ANISN/Pc generate the same fluxes for the same physical problem. This verifies 
that MGSEMI can be used as a tool for solving homogeneous semi-infinite media 
neutral particle transport problems. 
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4.2.3 Critical Slab Comparisons 

The FNCRIT and Anisn/Pc programs are used to determine critical slab 
thicknesses and resultant angular and scalar fluxes for two different cross sections 
and various discrete ordinates, Sn- The physical problem being analyzed is a single 
energy group, single finite slab whose thickness is determined by the amount of 
material necessary to achieve a self-sustaining neutron population. The material is 
defined by the total to scatter plus fission cross section ratio, c. 

Table 4.11 shows the critical slab thickness as calculated by the programs. 
Figures 4.6 through 4.9 show example angular and scalar comparison plots for the 
S32 and cross sections of 1.6 and 1.1. 

A normalization scheme is required to enable comparison between FNCRIT 
and Anisn/Pc. The normalization is performed in FNCRIT and sets its scalar flux 
at x = 0 equal to the ANISN/Pc scalar flux. The normalization parameter in 
FNCRIT is the ao parameter as discussed in Section A. 8.2. The analytical scalar 
flux at x = 0 for the Fn method is 

<M°) = 7^0 = 0anisn(O) (4. la) 

The scale parameter, a o, is solved for and used in the FNCRIT program. This 
requires the Anisn/Pc problem to be executed before the FNCRIT program. 

From the example plots shown and those not shown, FNCRIT and 
ANISN/Pc generate the same fluxes, for discrete ordinates of large order, for the 
same physical problem. This verifies that FNCRIT can be used as a tool for solving 
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critical thickness problems which is realized, will not be part of a spacecraft 
shielding problem. 

4.2.4 Scientific Literature Comparisons 

The last comparison is with calculations in the scientific literature. In 
Reference [12], two problems are analyzed with the Fn method. The first problem 
contains pathological cross sections. The cross sections contain a degenerate 
eigenvalue and a self scatter cross section of zero. A degenerate eigenvalue is a 
singularity problem encountered in the energy group particle transfer equations, 
equations (2.25), when s gg *v = i/q. The self scatter cross section being equal to 
zero creates another singularity in the energy group particle transfer equations, 
equations (2.24). The forerunner to MGSLAB is used to generate the fluxes for this 
cross section set. Because the degenerate eigenvalue derivation is dependent on a 
special case that will not manifest itself in most real material cross sections, the 
algorithm was not included in the current version of MGSLAB. The second 
problem is an analysis of gamma rays interacting with an iron slab. The original 
literature only included the angular fluxes at the slab boundaries; therefore, only 
these values are shown. Tables 4.12 through 4.15 show the pathological cross 
section boundary fluxes. Tables 4.16 through 4.19 show the iron slab boundary 
fluxes. These tables compare exactly, to a specified tolerance, with the fluxes in 
Reference [12]. This verifies that Fn method as implemented in this dissertation is 
proper and can be used to solve neutral particle transport problems like the 



galactic cosmic ray cascade. 



Table 4.1: Variation of Matrix Element Integral Quadrature to Determine the 

Number of Expansion Terms for Convergence for MGSLAB 


Integral 

Quadrature 

Expansion 

Terms 

Flux at n — 0.5 

x = 0.25 cm 

x = 0.50 cm 

x = 0.75 cm 

5 

51 

3.63639E-1 

5.70638E-1 

6.48887E-1 

10 

47 

3.63458E-1 

5.70738E-1 

6.49061E-1 

15 

47 

3.83459E-1 

5.70738E-1 

6.49060E-1 

25 

41 

3.63459E-1 

5. 70738 EM 

6.49060E-1 

50 

41 

3.63459E-1 

5.70738E-1 

6.49060E-1 

90 

41 

3.63459E-1 

5.70738E-1 

6.49060E-1 


Table 4.2: Variation of Matrix Element Integral Quadrature to Determine the 

Number of Expansion Terms for Convergence for MGSEMI 


Integral 

Quadrature 

Expansion 

Terms 

Flux at n = 0.5 

x = 0.25 cm 

x = 0.50 cm 

x = 0.75 cm 

5 

47 

1.29453 

2.26721 

2.98061 

10 

45 

1.29428 

2.26730 

2.98071 

15 

39 

1.29428 

2.26730 

2.98072 

25 

39 

1.29428 

2.26730 

2.98072 

50 

39 

1.29428 

2.26730 

2.98071 

90 

39 

1.29428 

2.26730 

2.98071 
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Table 4.3: Variation of Convergence Tolerance and the Resultant Fluxes for 

Mgslab 


Convergence 

Tolerance 

Expansion 

Terms 

Flux at /x = 0.5 

x = 0.25 cm 

x = 0.50 cm 

x = 0.75 cm 

1.0e-2 

13 

5.7864393E-02 

9. 0835148 E-02 

1.0329645E-01 

1.0e-3 

15 

5.7834357E-02 

9.0839556E-02 

1.0330132E-01 

1.0e-4 

23 

5.7847281E-02 

9.0836 134E-02 

1.0330116E-01 

1.0e-5 

31 

5.7846122E-02 

9. 0835786 E-02 

1.0330109E-01 

1.0e-6 

41 

5.7846242E-02 

9. 0835738 E-02 

1 .03301 10E-01 

1.0e-7 

59 

5.7846255E-02 

9.0835744E-02 

1 .03301 10E-0 1 

2.0e-8 

81 

5.7846256E-02 

9.0835744E-02 

1.0330110E-01 


Table 4.4: Variation of Convergence Tolerance and the Resultant Fluxes for 

Mgsemi 


Convergence 

Tolerance 

Expansion 

Terms 

Flux at /x = 0.5 

x = 0.25 cm 

x = 0.50 cm 

x = 0.75 cm 

1.0e-2 

9 

2.0583760E-01 

3.6091514E-01 

4.7444385E-01 

1.0e-3 

13 

2.0601616E-01 

3. 6085225 E-01 

4.7438955E-01 

1.0e-4 

21 

2.0598931E-01 

3. 6085700 E-01 

4.7439580E-01 

1.0e-5 

29 

2.0599167E-01 

3. 6085168 E-01 

4.7439558E-01 

1.0e-6 

39 

2.0599147E-01 

3.6085199 E-01 

4.7439553E-01 

1.0e-7 

49 

2.0599145E-01 

3. 6085198 E-01 

4.7439552E-01 

2.0e-8 

81 

2.0599145E-01 

3.6085198E-01 

4.7439552E-01 
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Table 4.5: 


Table 4.6: 


Variation of the Initial Number of Expansion Terms and the Number 
of Expansion Terms for Convergence for MGSLAB 


Initial Number 
of Expansion Terms 

Number of Expansion 
Terms for Convergence 

3 

29 

7 

29 

11 

29 

15 

29 

21 

29 


Variation of the Initial Number of Expansion Terms and the Number 
of Expansion Terms for Convergence for MGSEMI 


Initial Number 
of Expansion Terms 

Number of Expansion 
Terms for Convergence 

3 

29 

7 

29 

11 

29 

15 

29 

21 

29 
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Table 4.7: Variation of Basis Functions in Determining the Fluxes for MGSLAB 


P a (2f - 1) 
7 

Expansion 

Terms 

Flux at fi = 0.5 

x = 0.25 cm 

x = 0.50 cm 

x = 0.75 cm 

0.70 

33 

3.63459E-1 

5.70738E-1 

6.49060E-1 

0.75 

29 

3.63459E-1 

5.70738E-1 

6.49060E-1 

0.80 

31 

3.63459E-1 

5.70738E-1 

6.49060E-1 

0.85 

31 

3.63459E-1 

5.70738E-1 

6.49060E-1 

0.90 

31 

3.63459E-1 

5.70738E-1 

6.49060E-1 

0.95 

35 

3.63459E-1 

5.70738E-1 

6.49060E-1 

1.00 

41 

3.63459E-1 

5.70738E-1 

6.49060E-1 


Table 4.8: Variation of Basis Functions in Determining the Fluxes for MGSEMI 


P a { 2iF - 1) 

7 

Expansion 

Terms 

Flux at /i = 0.5 

x = 0.25 cm 

x — 0.50 cm 

x — 0.75 cm 

0.65 

33 

1.29428 

2.26730 

2.98071 

0.70 

27 

1.29428 

2.26730 

2.98071 

0.75 

29 

1.29428 

2.26730 

2.98071 

0.80 

27 

1.29428 

2.26730 

2.98071 

0.85 

27 

1.29428 

2.26730 

2.98071 

0.90 

31 

1.29428 

2.26730 

2.98071 

0.95 

33 

1.29428 

2.26730 

2.98071 

1.00 

39 

1.29428 

2.26730 

2.98071 





Table 4.9: Variation of the Number of Iterations on Interior Slab Interfaces and 

the Number of Expansion Terms for Convergence 


Number of Inner 
Iterations 

Expansion Terms 
for Convergence 

Flux at fi = 0.5 

x = 1.0 cm 

x = 2.0 cm 

x = 3.0 cm 

1 

31 

1.0247E-1 

1.5446E-2 

2.3006E-2 

2 

27 

1.0247E-1 

1.5446E-2 

2.3006E-2 

3 

25 

1.0247E-1 

1.5446E-2 

2.3006E-2 

4 

25 

1.0247E-1 

1.5446E-2 

2.3006E-2 

5 

23 

1.0247E-1 

1.5446E-2 

2.3006E-2 


Table 4.10: Comparison of Fluxes from a Four Slab Geometry versus a Single 

Slab with Three Interior Points at a Tolerance of 1.0 x 10 -6 


Geometry 

Flux at fi = 0.5 

x = 0.25cm 

x = 0.50cm 

x = 0.75cm 

i = 1.00cm 

4 Slabs 

3.63459E-1 

5.70738E-1 

6.49060E-1 

6.24782E-1 

1 Slab 

3.63459E-1 

5.70738E-1 

6.49060E-1 

6.24782E-1 





Table 4.11: 


Critical Slab Thickness for the FNCRIT and Anisn/Pc Programs 


Quadrature and Cross Section 

FNCRIT in cm 

Anisn/Pc in cm 

Sg and c = 1.6 

1.02392597657862 

1.04069 

Si 6 and c = 1.6 

1.02392597657862 

1.02701 

S 32 and c = 1.6 

1.02392597657862 

1.02476 

S 32 and c = 1.1 

4.22661933230336 

4.22834 
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Table 4.12: 


Left Boundary Angular Flux Values for the Pathological Cross 
Section Problem 


g 


0.05 

0.1 

0.2 

0.3 

0.4 

0.5 

1 

1.2813(-1) 1 

1.1660(-1) 

1 .0041 (- 1 ) 

8.8866(-2) 

7.9986(-2) 

7.2858(-2) 

2 

7.7363(-2) 

7.1154(-2) 

6.2123(-2) 

5.5490(-2) 

5.0281(-2) 

4.6037(-2) 

3 

5.6265(-2) 

5.2399(-2) 

4.6539(-2) 

4.2071(-2) 

3.8471(-2) 

3.5481(-2) 

4 

4.5318(-2) 

4.2624(-2) 

3.8381(-2) 

3.5035(-2) 

3.2276(-2) 

2.9946(-2) 

5 

3.8517(-2) 

3.6527(-2) 

3.3272(-2) 

3.0621 (-2) 

2.8389(-2) 

2.6473(-2) 

6 

3.3837(-2) 

3.2315(-2) 

2.9729(-2) 

2.7556(-2) 

2.5689(-2) 

2.4062(-2) 

7 

3.0396(-2) 

2.9207(-2) 

2.7105(-2) 

2.5284(-2) 

2.3687(-2) 

2.2276(-2) 

8 

2.7746(-2) 

2.6806(-2) 

2.5072(-2) 

2.3521(-2) 

2.2134(-2) 

2.0891(-2) 

9 

2.5634(-2) 

2.4886(-2) 

2.3442(-2) 

2.2106(-2) 

2.0887(-2) 

1.9779(-2) 

10 

2.3904(-2) 

2.3310(-2) 

2.2099(-2) 

2.0939(-2) 

1.9859(-2) 

1.8864(-2) 

11 

2.2458(-2) 

2.1988(-2) 

2.0970(-2) 

1.9957(-2) 

1.8994(-2) 

1.8094(-2) 

12 

2.1228(-2) 

2.0861(-2) 

2.0004(-2) 

1.9116(-2) 

1.8253(-2) 

1.7435(-2) 

13 

2.0166(-2) 

1.9885(-2) 

1.9166(-2) 

1.8385(-2) 

1.7609(-2) 

1.6862(-2) 

14 

1.9238(-2) 

1 .9031 (-2) 

1.8430(-2) 

1.7743(-2) 

1.7042(-2) 

1.6359(-2) 

15 

1.7895(-5) 

1.7895(-5) 

1.7895(-5) 

1.7895(-5) 

1.7894(-5) 

1.7894(-5) 

16 

9.1193(-3) 

9.1021(-3) 

8.9941(-3) 

8.8398(-3) 

8.6638(-3) 

8.4777(-3) 


1 Read as 1.2813 x 10 1 




Table 4.13: Left Boundary Angular Flux Values for the Pathological Cross 

Section Problem — Continued 


g 

V- 

0.6 

0.7 

0.8 

0.9 

1.0 

1 

6.6973(-2) 

6.2013(-2) 

5.7767(-2) 

5.4084(-2) 

5.0857(-2) 

2 

4.2494(-2) 

3.9480(-2) 

3.6881(-2) 

3.4613(-2) 

3.2615(-2) 

3 

3.2946(-2) 

3.0763(-2) 

2.8861(-2) 

2.7186(-2) 

2.5699(-2) 

4 

2.7943(-2) 

2.6199(-2) 

2.4666(-2) 

2.3306(-2) 

2.2090(-2) 

5 

2.4806(-2) 

2.3341 (-2) 

2.2042(-2) 

2.0881(-2) 

1.9838(-2) 

6 

2.2631(-2) 

2.1361(-2) 

2.0226(-2) 

1.9206(-2) 

1.8284(-2) 

7 

2. 1021 (-2) 

1.9897(-2) 

1.8886(-2) 

1.7971(-2) 

1.7140(-2) 

8 

1.9773(-2) 

1.8764(-2) 

1.7850(-2) 

1.7019(-2) 

1.6260(-2) 

9 

1.8773(-2) 

1.7858(-2) 

1.7023(-2) 

1.6259(-2) 

1.5558(-2) 

10 

1 .795 1 (-2) 

1.71 13(-2) 

1.6344(-2) 

1.5637(-2) 

1.4985(-2) 

11 

1.7260(-2) 

1.6488(-2) 

1.5775(-2) 

1.51 16(-2) 

1.4506(-2) 

12 

1.6669(-2) 

1.5954(-2) 

1.5290(-2) 

1.4673(-2) 

1.4100(-2) 

13 

1.6156(-2) 

1.5492(-2) 

1 .4871 (-2) 

1.4291 (-2) 

1.3749(-2) 

14 

1.5705(-2) 

1.5086(-2) 

1.4503(-2) 

1.3956(-2) 

1.3441(-2) 

15 

1.7894(-5) 

1.7894(-5) 

1.7893(-5) 

1.7893(-5) 

1.7893(-5) 

16 

8.2879(-3) 

8.0983(-3) 

7.9110(-3) 

7.7276(-3) 

7.5489(-3) 




Table 4.14: Right Boundary Angular Flux Values for the Pathological Cross 

Section Problem 


g 

V- 

0.05 

0.1 

0.2 

0.3 

0.4 

0.5 

1 

4.5701(-7) 

4.8998(-7) 

5.6471(-7) 

6.5894(-7) 

7.8516(-7) 

9.6666(-7) 

2 

7.3831 (-7) 

7.9749(-7) 

9.2983(-7) 

1.0958(-6) 

1.3185(-6) 

1.6409(-6) 

3 

9.6077(-7) 

1.0392(-6) 

1.2115(-6) 

1.4230(-6) 

1.7003(-6) 

2.0893(-6) 

4 

1.1809(-6) 

1.2788(-6) 

1.4910(-6) 

1.7477(-6) 

2.0786(-6) 

2.5330(-6) 

5 

1.4023(-6) 

1.5200(-6) 

1.7726(-6) 

2.0744(-6) 

2.4585(-6) 

2.9771(-6) 

6 

1.6263(-6) 

1.7641(-6) 

2.0575(-6) 

2.4046(-6) 

2.8416(-6) 

3.4235(-6) 

7 

1.8529(-6) 

2.01 12(-6) 

2.3460(-6) 

2.7385(-6) 

3.2281 (-6) 

3.8724(-6) 

8 

2.0822(-6) 

2.2612(-6) 

2.6377(-6) 

3.0759(-6) 

3.6177(-6) 

4.3236(-6) 

9 

2.3135(-6) 

2.5137 (-6) 

2.9322(-6) 

3.4160(-6) 

4.0098(-6) 

4.7762(-6) 

10 

2.5465(-6) 

2.7679(-6) 

3.2287(-6) 

3.7581(-6) 

4.4034(-6) 

5.2293(-6) 

11 

2.7806(-6) 

3.0234(-6) 

3.5265(-6) 

4.1013(-6) 

4.7976(-6) 

5.6820(-6) 

12 

3.0152(-6) 

3.2794(-6) 

3.8250(-6) 

4.4449(-6) 

5.1916(-6) 

6.1333(-6) 

13 

3.2497(-6) 

3.5354(-6) 

4.1232(-6) 

4.7880(-6) 

5.5843(-6) 

6.5821 (-6) 

14 

3.4835(-6) 

3.7907(-6) 

4.4206(-6) 

5.1297(-6) 

5.9749(-6) 

7.0275(-6) 

15 

2.3398(-9) 

2.3402(-9) 

2.3410(-9) 

2.3419(-9) 

2.3427(-9) 

2.3435(-9) 

16 

1.2655(-6) 

1.3423(-6) 

1.4865(-6) 

1.6302(-6) 

1.7799(-6) 

1.9399(-6) 
















































Table 4.15: 


Right Boundary Angular Flux Values for the Pathological Cross 
Section Problem — Continued 


g 

A* 

0.6 

0.7 

0.8 

0.9 

1.0 

1 

1.2644(-6) 

1.9188(-6) 

3.8488(-6) 

9.6935(-6) 

2.5162(-5) 

2 

2.1567(-6) 

3.0745(-6) 

4.7969(-6) 

7.9539(-6) 

1.3354(-5) 

3 

2.6827(-6) 

3.6719(-6) 

5.4072(-6) 

8.4198(-6) 

1.3377(-5) 

4 

3.2046(-6) 

4.2769(-6) 

6.0732(-6) 

9.0719(-6) 

1.3866(-5) 

5 

3.7257(-6) 

4.8834(-6) 

6.7563(-6) 

9.7887(-6) 

1.4525(-5) 

6 

4.2475(-6) 

5.4902(-6) 

7.4451(-6) 

1.0532(-5) 

1.5262(-5) 

7 

4.7704(-6) 

6.0965(-6) 

8.1350(-6) 

1.1287(-5) 

1.6036(-5) 

8 

5.2938(-6) 

6.7014(-6) 

8.8230(-6) 

1.2045(-5) 

1.6828(-5) 

9 

5.8170(-6) 

7.3039(-6) 

9.5070(-6) 

1.2799(-5) 

1.7625(-5) 

10 

6.3390(-6) 

7.9030(-6) 

1.0185(-5) 

1.3548(-5) 

1.8421(-5) 

11 

6.8589(-6) 

8.4974(-6) 

1.0857(-5) 

1.4288(-5) 

1.9209(-5) 

12 

7.3754(-6) 

9.0859(-6) 

1.1519(-5) 

1.5018(-5) 

1.9985(-5) 

13 

7.8876(-6) 

9.6676(-6) 

1.2172(-5) 

1.5734(-5) 

2.0749(-5) 

14 

8.3944(-6) 

1.0241 (-5) 

1.2813(-5) 

1.6437(-5) 

2.1495(-5) 

15 

2.3443(-9) 

2.3450(-9) 

2.3458(-9) 

2.3466(-9) 

2.3474(-9) 

16 

2.1142(-6) 

2.3070(-6) 

2.5235(-6) 

2.7700(-6) 

3.0545(-6) 




Table 4.16: Left Boundary Angular Flux Values for the Iron Slab Problem 


g 


0.05 

0.1 

0.2 

0.3 

0.4 

0.5 

1 

2.6676(-2) 

2.3995(-2) 

2.0336(-2) 

1.7801 (-2) 

1.5889(-2) 

1.4378(-2) 

2 

5.21 1 1 (-2) 

4.7223(-2) 

4.0437(-2) 

3.5657(-2) 

3.2006(-2) 

2.9093(-2) 

3 

3.7432(-2) 

3.4273(-2) 

2.9767(-2) 

2.6511 (-2) 

2.3980(-2) 

2.1931(-2) 

4 

3.6376(-2) 

3.3630(-2) 

2.9604(-2) 

2.6618(-2) 

2.4253(-2) 

2.2312(-2) 

5 

3.5820(-2) 

3.3465(-2) 

2.9892(-2) 

2.7159(-2) 

2.4947(-2) 

2.3100(-2) 

6 

3.6063(-2) 

3.4070(-2) 

3.0917(-2) 

2.841 1(-2) 

2.6329(-2) 

2.4556(-2) 

7 

2.4795(-2) 

2.3650(-2) 

2.1757(-2) 

2.0194(-2) 

1.8861(-2) 

1.7704(-2) 

8 

2.6132(-2) 

2.5122(-2) 

2.3376(-2) 

2.1882(-2) 

2.0576(-2) 

1.9421(-2) 

9 

2.8500(-2) 

2.7614(-2) 

2.5995(-2) 

2.4547(-2) 

2.3245(-2) 

2.2069(-2) 

10 

3.2597(-2) 

3.1830(-2) 

3.0318(-2) 

2.8889(-2) 

2.7560(-2) 

2.6330(-2) 

11 

3.9796(-2) 

3.9167 (-2) 

3.7765(-2) 

3.6333(-2) 

3.4940(-2) 

3.3610(-2) 

12 

5.1910(-2) 

5.1584(-2) 

5.0494(-2) 

4.9168(-2) 

4.7762(-2) 

4.6346(-2) 

13 

1.5781(-2) 

1.6779(-2) 

1.7971(-2) 

1.8609(-2) 

1.8935(-2) 

1.9064(-2) 

14 

1.3336(-2) 

1.4244(-2) 

1.5458(-2) 

1.6229(-2) 

1.6731(-2) 

1.7049(-2) 

15 

8.8737(-3) 

9.4910(-3) 

1.0406(-2) 

1.1073(-2) 

1.1575(-2) 

1.1959(-2) 

16 

2.3930(-3) 

2.5592(-3) 

2.8282(-3) 

3.0458(-3) 

3.2278(-3) 

3.3821(-3) 

17 

6.0379(-4) 

6.4190(-4) 

7.0559(-4) 

7.5958(-4) 

OO 

o 

o 

0 

1 

8.4932(-4) 

18 

4.0769(-5) 

4.3372(-5) 

4.7756(-5) 

5.1532(-5) 

5.4913(-5) 

5.7999(-5) 

19 

6.1145(-6) 

6.4725(-6) 

7.0782(-6) 

7.6032(-6) 

8.0771 (-6) 

8.5134(-6) 














































































Table 4.17: Left Boundary Angulax Flux Values for the Iron Slab Problem 

Continued 


g 


0.6 

0.7 

0.8 

0.9 

1.0 

1 

1.3145(-2) 

1.2117 (-2) 

1.1244(-2) 

1.0492(-2) 

9.8369(-3) 

2 

2.6698(-2) 

2.4687(-2) 

2.2970(-2) 

2.1485(-2) 

2.0185(-2) 

3 

2.0228(-2) 

1.8784(-2) 

1.7541(-2) 

1.6459(-2) 

1.5507(-2) 

4 

2.0679(-2) 

1.9282(-2) 

1.8071(-2) 

1.7008(-2) 

1.6067(-2) 

5 

2.1526(-2) 

2.0165(-2) 

1.8972(-2) 

1.7918(-2) 

1.6979(-2) 

6 

2.3021(-2) 

2.1676(-2) 

2.0485(-2) 

1.9422(-2) 

1.8467(-2) 

7 

1.6688(-2) 

1.5785(-2) 

1.4978(-2) 

1.4251(-2) 

1.3592(-2) 

8 

1.8392(-2) 

1.7468(-2) 

1.6633(-2) 

1.5874(-2) 

1.5183(-2) 

9 

2.1004(-2) 

2.0035(-2) 

1.9150(-2) 

1.8338(-2) 

1.7592(-2) 

10 

2.5194(-2) 

2.4144(-2) 

2.3173(-2) 

2.2273(-2) 

2.1437(-2) 

11 

3.2352(-2) 

3.1167(-2) 

3.0054(-2) 

2.9009(-2) 

2.8027(-2) 

12 

4.4954(-2) 

4.3603(-2) 

4.2303(-2) 

4.1057(-2) 

3.9867(-2) 

13 

1.9060(-2) 

1.8967(-2) 

1 .88 10(-2) 

1.8608(-2) 

1.8376(-2) 

14 

1.7237(-2) 

1.7329(-2) 

1.7348(-2) 

1.7312(-2) 

1.7235(-2) 

15 

1.2252(-2) 

1.2473(-2) 

1.2636(-2) 

1.2753(-2) 

1.2832(-2) 

16 

3.5139(-3) 

3.6268(-3) 

3.7236(-3) 

3.8067(-3) 

3.8779(-3) 

17 

8.8743(-4) 

9.2194(-4) 

9.5330(-4) 

9.8188(-4) 

1.0080(-3) 

18 

6.0844(-5) 

6.3487(-5) 

6.5953(-5) 

6.8263(-5) 

7.0432(-5) 

19 

8.9197(-6) 

9.3012(-6) 

9.6613(-6) 

1.0003(-5) 

1.0327(-5) 



Table 4.18: Right Boundary Angular Flux Values for the Iron Slab Problem 




0.05 

0.1 

0.2 

0.3 

0.4 


1 

4.6286(-5) 

4.9431(-5) 

5.7083(-5) 

6.7649(-5) 

8.9146(-5) 

1.7505(-4) 

2 

9.8793(-5) 

1.0537(-4) 

1.2091(-4) 

1.4138(-4) 

1.7025(-4) 

2.1330(-4) 

3 

8.3317(-5) 

8.8923(-5) 

1.0170(-4) 

1.1787(-4) 

1.3969(-4) 

1.7073(-4) 

4 

9.3154(-5) 

9.9491(-5) 

1.1351(-4) 

1.3070(-4) 

1.5305(-4) 

1.8364(-4) 

5 

1.0571(-4) 

1.1298(-4) 

1.2863(-4) 

1.4722(-4) 

1.7059(-4) 

2.0138(-4) 

6 

1.2265(-4) 

1.3119(-4) 

1.4907(-4) 

1.6967(-4) 

1.9472(-4) 

2.2656(-4) 

7 

9.4413(-5) 

1.0103(-4) 

1.1460(-4) 

1.2983(-4) 

1.4790(-4) 

1.7019(-4) 

8 

1.0859(-4) 

1.1624(-4) 

1.3167(-4) 

1.4865(-4) 

1.6836(-4) 

1.9212(-4) 

9 

1.2831 (-4) 

1.3739(-4) 

1.5538(-4) 

1.7478(-4) 

1.9680(-4) 

2.2274(-4) 

10 

1.5735(-4) 

1.6847(-4) 

1.9015(-4) 

2.1301(-4) 

2.3840(-4) 

2.6758(-4) 

11 

2.0326(-4) 

2.1750(-4) 

2.4478(-4) 

2.7287(-4) 

3.0332(-4) 

3.3740(-4) 

12 

2.81 11 (-4) 

3.0041(-4) 

3.3661(-4) 

3.7286(-4) 

4.1 102(-4) 

4.5244(-4) 

13 

1.6780(-4) 

1.8023(-4) 

2.0290(-4) 

2.2480(-4) 

2.4706(-4) 

2.7036(-4) 

14 

1.5666(-4) 

1.6825(-4) 

1.8907(-4) 

2.0873(-4) 

2.2820(-4) 

2.4801(-4) 

15 

1.1581(-4) 

1.2400(-4) 

1.3847(-4) 

1.5181(-4) 

1.6465(-4) 

1.7733(-4) 

16 

3.6635(-5) 

3.9089(-5) 

4.3366(-5) 

4.7232(-5) 

5.0873(-5) 

5.4381 (-5) 

17 

9.5847(-6) 

1.0174(-5) 

1.1189(-5) 

1.2093(-5) 

1.2930(-5) 

1.3722(-5) 

18 

6.6624(-7) 

7.0823(-7) 

7.7968(-7) 

8.4231 (-7) 

8.9961(-7) 

9.5318(-7) 

19 

1.0024(-7) 

1.0606(-7) 

1.1596(-7) 

1.2461(-7) 

1.3249(-7) 

1.3983(-7) 































































Table 4.19: Right Boundary Angular Flux Values for the Iron Slab Problem 

Continued 


g 


0.6 

0.7 

0.8 

0.9 

1.0 

1 

4.8010(-4) 

1.2490(-3) 

2.7431(-3) 

5.1717(-3) 

8.6617 (-3) 

2 

2.7680(-4) 

3.6438(-4) 

4.7562(-4) 

6.0693(-4) 

7.5306(-4) 

3 

2.1516(-4) 

2.7581(-4) 

3.5305(-4) 

4.4501(-4) 

5.4851 (-4) 

4 

2.2618(-4) 

2.8350(-4) 

3.5647(-4) 

4.4389(-4) 

5.4319(-4) 

5 

2.4288(-4) 

2.9791(-4) 

3.6772(-4) 

4.5175(-4) 

5.4802(-4) 

6 

2.6810(-4) 

3.2208(-4) 

3.9012(-4) 

4.7223(-4) 

5.6704(-4) 

7 

1.9848(-4) 

2.3455(-4) 

2.7964(-4) 

3.3408(-4) 

3.9731(-4) 

8 

2.2158(-4) 

2.5850(-4) 

3.0422(-4) 

3.5935(-4) 

4.2364(-4) 

9 

2.5414(-4) 

2.9272(-4) 

3.3996(-4) 

3.9671(-4) 

4.6309(-4) 

10 

3.0202(-4) 

3.4340(-4) 

3.9329(-4) 

4.5287(-4) 

5.2261 (-4) 

11 

3.7652(-4) 

4.2228(-4) 

4.7635(-4) 

5.4016(-4) 

6.1467(-4) 

12 

4.9844(-4) 

5.5048(-4) 

6.1015(-4) 

6.7904(-4) 

7.5855(-4) 

13 

2.9526(-4) 

3.2232(-4) 

3.5212(-4) 

3.8526(-4) 

4.2228(-4) 

14 

2.6857(-4) 

2.9022(-4) 

3.1330(-4) 

3.3815(-4) 

3.6512(-4) 

15 

1.9006(-4) 

2.0300(-4) 

2.1628(-4) 

2.3004(-4) 

2.4440(-4) 

16 

5.7808(-5) 

6.1193(-5) 

6.4563(-5) 

6.7942(-5) 

7.1348(-5) 

17 

1.4483(-5) 

1.5219(-5) 

1.5937(-5) 

1.6641(-5) 

1.7335(-5) 

18 

1.0039(-6) 

1.0525(-6) 

1.0992(-6) 

1.1445(-6) 

1.1885(-6) 

19 

1.4675(-7) 

1.5333(-7) 

1.5965(-7) 

1.6573(-7) 

1.7161(-7) 
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Direction At 

Figure 4.1: ANISN/Pc versus MGSLAB Angular Fluxes for Three Regions and 

Three Energy Groups with a Beam Source at Slab Positions of 
0.0 cm, 0.1 cm, 0.2 cm, 0.4 cm, 0.6 cm, 0.8 cm, and 1.0 cm 



Figure 4.2: ANISN/Pc versus MGSLAB Scalar Fluxes for Three Regions and 

Three Energy Groups with a Beam Source 
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Direction $jl 


Figure 4.3: ANISN/Pc versus MGSLAB Angular Fluxes for Three Regions and 

Three Energy Groups with a Cubic Spline Fit Distributed Source at 
Slab Positions of 0.0 cm, 0.1 cm, 0.2 cm, 0.4 cm, 0.6 cm, 0.8 cm, and 
1.0 cm 



Figure 4.4: ANISN/Pc versus MGSLAB Scalar Fluxes for Three Regions and 

Three Energy Groups with a Cubic Spline Fit Distributed Source 
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Figure 4.5: Anisn/Pc versus MGSEMI Boundary Angular Fluxes for Three 

Energy Group with a Beam Source at the Boundary 






Scalar Flux 
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Figure 4.6: ANISN/Pc versus FNCRIT Angular Flux Values for an S32 Problem 

and c = 1.6 



Figure 4.7: Anisn/Pc versus FNCRIT Scalar Flux Values for an S32 Problem and 

c = 1.6 
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Chapter 5 
Results 


In Chapter 4, it was shown that the MGSLAB, MGSEMI, and FNCRIT 
programs generate verifiable results for beam, isotropic, distributed source, and 
critical slab problems. The MGSEMI and MGSLAB programs are coupled to the 
galactic ion transport program, GlT, and applied to the galactic cosmic ray 
cascade problem. The specific cross section, multiplicity, path length, and stopping 
power models used in the coupled programs are discussed. For all calculations, the 
target material is aluminum which is composed of 100% Al-27 with a density of 


2.696 


5.1 Cross Section and Multiplicity Models 
Two cross section models are used in the galactic cosmic ray cascade 
problem: a nuclear liquid drop model for the GlT program and the Endf/B V 
cross sections for the MGSEMI and MGSLAB programs. The liquid drop model can 
be used in the GlT program because the speed of the interacting ions is assumed 


to be large. 
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5.1.1 GlT Cross Sections 

In the nuclear liquid drop model, it is assumed that the total nucleax charge 

is spread uniformly throughout the entire spherical nucleax volume and the nuclear 

density is constant. From these assumptions, the nucleax volume is proportional to 

the atomic mass, or V oc A. Therefore, the nucleax radius is proportional to the 

cube root of the atomic mass, or R = Rq x A 3 . The speed of the ions is assumed 

large; therefore, the cross sectional area of the ion nucleus is the cross section of 

interaction because resonance and quantum effects are negligible. The cross 

sectional area of the target nucleus is assumed identical for all target nuclei and is 

not included explicitly in the interaction cross section, but is included in the 

proportionality constant, Rq. For the liquid drop model, the cross sectional area is 

2 

proportional to the square of the radius, or crj oc A*. In the GlT model, it is 
assumed that all ions contain only protons. Thus, A is defined as the charge 
number, j. The cross section can then be written as 

(5.1) 

where, a is the proportionality constant or the cross section normalization 

parameter. From Reference [6], the value of a representative of an air shield is 
2 

0.01247 This value accounts for the cross sectional area of the target nuclei 
and the proportionality constant, Rq, specified above. 

To use these cross sections in the MGSEMI and MGSLAB programs, the 
values are multiplied by the density of the target material. In addition, a must be 
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computed for the target material. The cross section normalization parameter is 


proportional to the cross sectional axea of the target material, A 3 ^^. The 


normalization parameter of the new target material can be related to that of an air 


shield 


(5.2) 


^target — ^air 


Uarget 


L A 


air J 


Therefore, the value for <7 ai is 


a M = 0.01247 


27 


L 14.4 J 


3 „ , 

= 0.1896 


The above constants are used in the GlT, MGSEMI, and MGSLAB programs 
to simulate the galactic cosmic ray cascade. 


5.1.2 Git Multiplicities 

The GlT cross section model represents the ion interaction probability for 
the target nuclei. The multiplicity, describes the number of j th ions coming 
from each fragmentation of the A: th ion. The charge number of the k th ion must be 
conserved after the interaction; therefore 


Mm = 


k > j 


k - 1 
0 k < j. 


(5.3) 


If the resultant ion charge number is multiplied by the multiplicity term and 
summed for all resultant ions, the outcome is equal to the charge number of the 
incoming ion; that is, charge is conserved. 
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As previously noted, the GlT model does not account for neutrons; 
however, ail ions except hydrogen contain neutrons. The GlT model is modified 
through the multiplicity term to accept neutrons. This is accomplished by defining 
a fraction, /„, of ail particles created with a charge of one as neutrons. These 
neutrons are then given a charge of zero. The value of /„ for the analyses in this 
dissertation is set at 0.5. 

The new multiplicity term for the creation of protons from the 
fragmentation of ion k is 

M U = (Wnljlf. (5.3b) 

For neutrons, the multiplicity term is 

Mo,* = /nT~— r- (5.4a) 

k — 1 

The importance of the f n fraction will be discussed in the conclusions; however, it 
allows the original GlT model to be used as a starting point for the source used in 
the neutral particle transport model. 

5.1.3 MGSEMI and MGSLAB Cross Sections 

A simple and accurate analytical model does not exist for neutron cross 
sections in Al-27. However, extensive experimental data do exist up to neutron 
energies of 20 MeV and are provided in the ENDF/B V database. This database is 
used to generate eleven energy group neutron cross sections for Al-27 with down 
scattering only. The NJOY program (Reference [26]) accessed and collapsed the 
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database to create the desired cross sections. Initial ion beams of 20 MeV per 
nucleon are not typical for the galactic cosmic ray cascade. A typical value is 
1000 MeV per nucleon. In order to create a cross section set to match the 
problem, the Endf/B V energy group of 10 MeV to 20 MeV is extended to 
100 MeV. To create a twelfth energy group, 100 MeV to 1000 MeV, the 20 MeV 
cross sections were decreased slightly and used. For high energies, these extended 
cross sections do not represent anything physical, but the point here is to 
demonstrate behavior, not generate accurate numbers. Table 5.1 shows the energy 
groups and macroscopic cross sections generated. The eleven group cross section 
will be referred to as the limited cross section set and the twelve group cross 
section set will be referred to a s the extend cross section set. 



The path length and proton stopping power model used in the galactic 
cosmic ray cascade is based on a parametric form of the proton range, Rp(E), 
(Reference [27]) 

Rp(E) = a' Q ln (l + a[E n °) , (5.5a) 

where, the parametric constants are 

a' 0 = 500, = 3.66 x 10" 6 , and n' 0 = 1.79. (5.5b) 

Therefore, from the definition of path length 


s 


= s( E 0 ,E) = Rp(E 0 ) - Rp(E), 


(5.6) 



Ill 


where, Eo is the energy of a proton in the initial beam of ions. From this path 
length, a stopping power can be determined (called the Wilson stopping power), 


S?(E) = - 


dE 

ds 


1 + tt , 1 E n ° 

a' 0 a[n' 0 E"° 


ni-l ' 


(5.7) 


For ea.se in computation, the Wilson stopping power can be simplified to a linear 
function of energy by setting n' Q to one and changing the two parametric constants, 
ol 0 and Ofj, to reflect the original function value. 

The form of this new simplified stopping power, Sp, is 


Sp(E) = <zE + 6 (5.8a) 

where, a and b are to be determined. Because the simplified stopping power is 
linear, the slope, a, can be extracted from two values of the stopping power known 
at two different energies. The two stopping powers are determined using the more 
accurate Wilson stopping power in equation (5.7). Therefore, the slope is 
Sp(Ei) — Sp(E 2 ) 

a = — — 

Ei — E2 

A simplified proton range function can be specified for the condition of 
linearity using the form of equation (5.5a) 



Rp(E) = a 0 ln(l + a x E) , 


(5.9) 


where, ao and o x are the new parametric constants with values that must reflect 


the original function. 
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Using this new proton range, the simplified path length is 

s = Rp(E 0 ) — a 0 ln(l + aiE) . (5.10) 

The simplified stopping power is determined from the simplified path length and 
compared with equation (5.8a) 

Sp(E) = ~ = — E + — ; (5.11a) 

as QLq OoGfi 

therefore, 

a 0 — - and b — . (5.11b) 

a <^o a i 

To determine the simplified proton range function is evaluated at Eo 

Rp(E 0 ) = a 0 ln(l + aiE 0 ) . (5.12a) 

If the proton range at Eo is evaluated using the more accurate formulation in 
equation (5.5a), then the above equation can be solved for a\ 

i _ pRp(Eo)/? 

a. = f . (5.12b) 

To extend this model so it can be used in the MGSEMI and MGSLAB 
programs, the energy and path length variables must be discretized as explained in 
Section 5.1.3. This only affects the nomenclature used in the above equations. The 
variable changes are 

E => Eg and s => s g . 

The simplified stopping power and path length models described above are 
used in the GlT, MGSEMI, and MGSLAB programs. 
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5.3 Application to the Galactic Cosmic Ray Cascade 
The galactic cosmic ray cascade algorithm will be demonstrated in two 
ways. The first will be to show examples of the algorithm. A beam of fluorine ions 
incident on a finite and semi-infinite aluminum slab will be analyzed. The second 
demonstration will show a study of the transmitted scalar flux for varying 
thickness of aluminum slabs. 

The specific format of the results are the neutron angular, and 

scalar, (f> g (x ), flux profiles, ion flux profiles, xf>j(x,E g ), and GlT source values, 
for each energy group. 

The neutron flux profiles, both angular and scalar, are the primary results 
for this work. From these profiles, new low energy neutron models can be added to 
BRYNTRN and verified. Once a suitable model exists for BRYNTRN, the neutron 
dose rates can then be calculated. The GlT, MGSEMI, and MGSLAB programs 
have been verified separately (Reference [6] and Chapter 4); therefore, these results 
are within a specified numerical tolerance within the given physical assumptions. 

5.3.1 Demonstration of the Coupled Algorithms 

To demonstrate the coupled programs, analyses of a 1000 MeV per nucleon 
fluorine beam incident on semi-infinite and finite aluminum slabs are performed. 

To determine the importance of the neutrons with energies above 20 MeV, the 
extended cross section set is replaced with the limited cross section set and the 
same analyses are performed again. 
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5. 3. 1.1 Semi-infinite and Finite Slab Analysis for the Extended 
Cross Section Set. A beam of fluorine ions at an energy of 1000 MeV per 
nucleon axe normally incident on a semi-infinite aluminum slab. The resultant ion 
flux profiles for energies of 999 MeV per nucleon, 100 MeV per nucleon, and 
10 MeV per nucleon are shown in Figures 5.1 through 5.3. There is no change in 
the energy group profiles below 10 MeV per nucleon, therefore they are not shown. 
The reason they are identical is that the path length or penetration depth, 
s(Eo,E), is logarithmic with respect to energy and the small energy values, 
compared to 1000 MeV, do not alter the path length enough to change the flux 
profiles. The next set of profiles, Figures 5.4 through 5.6, show the neutron source 
used in the Fn algorithm. 

For this example, the Fn algorithm’s convergence tolerance is 5.0 x 10~ 3 to 
allow group six to converge. To compensate for the low convergence tolerance, the 
initial number of expansion terms is set to forty-one so that all the other groups 
converged with a maximum relative error approximately 1.0 x 10“ 6 . Only two 
points are slow to converge in group six. This is deemed unimportant to the 
overall problem allowing the convergence tolerance and initial number of expansion 
terms fix mentioned above. The reason for the slow convergence rate is due to the 
formulation of the integral in equation (2.58b). When the interior edit point nears 
the integral value is small and the resultant source is not large enough to let 
the Fn algorithm converge. The basis functions are modified Legendre polynomials 
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with 7 equal to 0.95. The number of direction edit points is eight which are set to 
the zeros of the Legendre polynomials. With these parameters, the angular 
neutron flux profiles at the slab boundary, 100 centimeters, and 200 centimeters 
are shown in Figures 5.7 through 5.10. 

Figure 5.11 shows a three dimensional surface profile of the scalar flux by 
energy group and position. Forty, evenly spaced internal positions to a depth of 
200 centimeters allow a fine enough grid to see the behavior of the neutrons well 
into the media. 

The same set of profiles shown for the semi-infinite media are also shown 
for a finite slab. Figures 5.12 through 5.14 are the ion flux profiles from 1000 Mev, 
100 MeV, and 10 MeV per nucleon. Figures 5.15 through 5.17 show the neutron 
source used in the Fn algorithm. Some of the parameters are changed for the finite 
slab case. The slab thickness is 120 centimeters and only twenty internal slab 
positions are used to determine the scalar flux. The basis functions are the shifted 
Legendre polynomials. The convergence tolerance is lowered to 2.0 x 10 -4 and the 
initial number of expansion terms is lowered to thirty-one. From these parameters, 
the angular neutron flux profiles are shown in Figures 5.18 through 5.21. The 
three dimensional surface profile of the scalar flux is shown in Figure 5.22. 

5. 3. 1.2 Semi-infinite and Finite Slab Analysis for the Limited Cross Section Set. 
To determine if the created cross sections for energy group twelve are important, 
the same problems as described above are executed for the limited cross section set 




116 


created by NJOY. The first energy group now only represents neutrons with 
energies between 10 MeV and 20 MeV. The results for the semi-infinite media are 
shown in Figures 5.23 through 5.33. The results for the finite slab are shown in 
Figures 5.34 through 5.44. 

5. 3. 1.3 Discussion. This section describes and discusses the profiles identified 
above. Also, the differences between the limited and extended cross section sets 
are discussed. 

A detailed discussion of the ion flux profiles in Figures 5.1 through 5.3, 
Figures 5.12 through 5.14, Figures 5.23 through 5.25, and Figures 5.34 
through 5.36, is given in Reference [6]; however, several important features of the 
profiles are highlighted here. Based on the continuous slowing down 
approximation, as an ion in the incident beam loses a given amount of energy, 

(E 0 — Ei), the ion traverses a known path length, s(E 0 , Ei) or s\. If the ion has 
not collided with a target nuclei, the ion travels a distance of xj = The 
number of ions from the initial beam that have not collided with the target and 
have slowed down to an energy of Ei is e~ (TjXj . The number of Ei ions that have 
fragmented into lighter ions (0 < j < J - 1) at xj is crj'&jix j , Et) or a J e ~ <7jXj . 
These secondary ions travel a maximum distance of Xj = at an energy of Ej. 
The source of J — 1 ions can only begin to generate these ions at xj. Because 
there is only a point source for J — 1 ions, the flux profile at energy Ei decays in 
the range xj to xj_ x . In this region, the J — 1 ions continuously generate lighter 
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ions. As the J — 1 ion flux decays, the number of J — 2 ions increases. However, 
when the source of the J — 2 ions abates, the J — 2 ion flux profile decays like the 
J — 1 ion profile until it reaches its maximum range. At xj_ i, the flux profile for 
J — 2 ions has an infinite slope change. This indicates that the source of J — 2 ions 
(i.e., those produced by J — 1 ions only) has gone abruptly to zero. This 
discontinuity occurs in all flux profiles, but the flux profiles become smoother 1 for 
the lighter ions because the discontinuous source at ion J (the original cause of the 
discontinuity) is a smaller fraction of the total source of the lighter ions. These 
features can be seen in varying detail in all of the ion flux profile plots. 

The neutron source profiles in Figures 5.4 through 5.6, Figures 5.15 
through 5.17, Figures 5.26 through 5.28, and Figures 5.37 through 5.39, have the 
same discontinuity as the ion flux profiles except the discontinuities are 
superimposed on one another throughout the slab because the neutron source at a 
certain energy is the summation of all higher energy ion sources. If the leading 
edge of the neutron source profile from Figure 5.4 is magnified, as in Figure 5.45, 
then the discontinuities can be seen in detail. The discontinuities still occur at -4*-, 
but the discontinuities for all ions at energy levels above and including group g are 
superimposed upon one another. The smoothing of the source profiles, as seen in 
the ion profiles, is also a feature of the neutron source profiles, and is enhanced for 
the larger energy groups. This enhancement is a result of adding the source from 


x The discontinuity manifests itself in higher and higher derivatives of the ion flux. 
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all ions and all higher energy groups together. In the regions of the slab where ions 
are generating neutrons, the number of neutrons increases. As the source of 
neutrons stops, the neutrons then start to decrease due to scattering and 
absorption. 

The neutron angular flux profiles in Figures 5.7 through 5.10, Figures 5.18 
through 5.21, Figures 5.29 through 5.32, and Figures 5.40 through 5.43, show that 
the angular fluxes are fairly constant with angle in the slab interiors. If there is a 
direction preference, then it is in the forward direction, the original direction of the 
ion beam. The transmitted boundary fluxes for the finite slabs show this forward 
peak preference to an even greater extent than the interior fluxes. The boundary 
fluxes for the semi-infinite slabs and the left face reflective fluxes for the finite 
slabs show a backward peaked flux because of the relative size difference between 
the aluminum nucleus and the neutron. For the finite slabs, the transmitted fluxes 
within each energy group are approximately the same value. This is due to the 
isotropic scattering assumption. 

The neutron scalar flux profiles in Figures 5.11, 5.22, 5.33, and 5.44 show 
that the scalar neutron flux peaks at about fifty centimeters into the slabs in the 
energy group that represents 0.1 MeV. The placement of the flux peaks within the 
slab matches that of the source; however, the source is largest in group three. 

There is a difference between the fluxes for semi-infinite and finite slabs. 
The angular fluxes are generally larger within the finite slabs which is due to the 
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size of the semi-infinite slab. The neutrons have more region to fill, so the numbers 
of a neutrons at each angle and position is smaller than if a boundary is present. 

Upon comparison of the limited cross section fluxes to the extended cross 
section fluxes (ignoring the first group), the limited fluxes are smaller by a factor of 
two in the interior and a factor of ten at the boundaries. This emphasizes the need 
for higher energy cross sections. This will discussed in detail in the conclusions. 

5.3.2 Slab Thickness Study 

As a practical example in using this algorithm, a study of aluminum slab 
thickness versus transmitted scalar flux, which can be related to exposure or dose, 
is undertaken. A convergence tolerance of 1G -4 is used for the extended cross 
section set to generate scalar fluxes at the transmission boundary for slabs ranging 
from thirty-five centimeters to 150 centimeters. Figure 5.46 shows the transmitted 
fluxes versus slab thickness. It shows an atypical peak around eighty centimeters. 
Common sense would dictate that as the slab gets thicker, the flux would decrease; 
however, the physics of the galactic cosmic ray cascade dictates that neutrons are 
created within the slab instead of at the slab boundary. If the slab is thin, then 
not enough material exists to create neutrons. If the slab is thick, then enough 
material exists to shield the neutrons created inside it. If the slab is as thick as 
where the neutron source peaks, then the transmitted neutron flux is at its 
highest. The implications to spacecraft design are very important. Certain wall 
thicknesses could produce larger neutron doses that thinner or thicker walls. The 
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thickness and material of the shield wall is a critical component to occupant safety 
which is the reason for the interest in this topic. 
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Table 5.1: Twelve Group, Extended, Down Scatter Al-27 Macroscopic Neutron 

Cross Sections Generated from Endf/B V and NJOY 


g Energy Total 

MeV o a 



msssmmmsM 

4.40101-2 

1E+2 

1.05095-1 

4.06120-2 


1.74551-1 

8.35890-3 

1E+0 

2.47382-1 

6.14990-4 

IE-1 

2.80333-1 

6.85712-6 

IE-2 

8.97516-2 

9.72033-8 

IE 3 

8.13180-2 

1.97605-9 

IE-4 

8.15770-2 

0.0 

IE-5 

8.24746-2 

0.0 

IE-6 

8.61492-2 

0.0 

IE-7 

9.23780-2 

0.0 

IE-8 

1.17601-1 

0.0 
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Figure 5.1: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 999.0 MeV (Energy Group 1) for a Semi-infinite 
Aluminum Media using the Extended Cross Section Set 



Figure 5.2: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 100.0 MeV (Energy Group 2) for a Semi-infinite 
Aluminum Media using the Extended Cross Section Set 
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Figure 5.3: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 10.0 MeV (Energy Group 3) for a Semi-infinite 
Aluminum Media using the Extended Cross Section Set 



Position in cm 


Figure 5.4: GlT Neutron Source Profile for Energy Groups One through Four for 

a Fluorine Beam of 1000 MeV per Nucleon and a Semi-infinite 
Aluminum Media using the Extended Cross Section Set 
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Figure 5.5: GlT Neutron Source Profile for Energy Groups Five through Eight 

for a Fluorine Beam of 1000 MeV per Nucleon and a Semi-infinite 
Aluminum Media using the Extended Cross Section Set 



Figure 5.6: GlT Neutron Source Profile for Energy Groups Nine through Twelve 

for a Fluorine Beam of 1000 MeV per Nucleon and a Semi-infinite 
Aluminum Media using the Extended Cross Section Set 
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Figure 5.7: Angular Neutron Flux Profile for x = 0 cm, x = 100 cm, and 

x = 200 cm for Energy Groups One through Three, a Fluorine Beam 
of 1000 MeV per Nucleon, and a Semi-infinite Aluminum Media using 
the Extended Cross Section Set 
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Figure 5.8: Angular Neutron Flux Profile for x = 0 cm, x = 100 cm, and 

x = 200 cm for Energy Groups Four through Six, a Fluorine Beam of 
1000 MeV per Nucleon, and a Semi-infinite Aluminum Media using 
the Extended Cross Section Set 
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Figure 5.9: Angular Neutron Flux Profile for x = 0 cm, x = 100 cm, and 

x = 200 cm for Energy Groups Seven through Nine, a Fluorine Beam 
of 1000 MeV per Nucleon, and a Semi-infinite Aluminum Media using 
the Extended Cross Section Set 







Neutron Directional Flux 
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Figure 5.11: Scalar Neutron Flux Profile by Energy Group with a Fluorine Beam 

of 1000 MeV per Nucleon and a Semi-infinite Aluminum Media 
using the Extended Cross Section Set 
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Figure 5.12: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 999.0 MeV (Energy Group 1) for a Finite 
Aluminum Slab using the Extended Cross Section Set 



Figure 5.13: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 100.0 MeV (Energy Group 2) for a Finite 
Aluminum Slab using the Extended Cross Section Set 
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Figure 5.14: Git Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 10.0 MeV (Energy Group 3) for a Finite Aluminum 
Slab using the Extended Cross Section Set 



Figure 5.15: Git Neutron Source Profile for Energy Groups One through Four 

for a Fluorine Beam of 1000.0 MeV per Nucleon and a Finite 
Aluminum Slab using the Extended Cross Section Set 
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Figure 5.16: GlT Neutron Source Profile for Energy Groups Five through Eight 

for a Fluorine Beam of 1000 MeV per Nucleon and a Finite 
Aluminum Slab using the Extended Cross Section Set 



Figure 5.17: GlT Neutron Source Profile for Energy Groups Nine through Twelve 

for a Fluorine Beam of 1000 MeV per Nucleon and a Finite 
Aluminum Slab using the Extended Cross Section Set 








Direction /i 


Figure 5.19: Angular Neutron Flux Profile for x = 0 cm, x = 60 cm, and 

x = 120 cm for Energy Groups Four through Six, a Fluorine Beam 
of 1000 MeV per Nucleon, and a Finite Aluminum Slab using the 
Extended Cross Section Set 





135 



0.00 0.25 0.50 0.75 1.00 

Direction 11 


Figure 5.20: Angular Neutron Flux Profile for x = 0 cm, x = 60 cm, and 

x = 120 cm for Energy Groups Seven through Nine, a Fluorine 
Beam of 1000 MeV per Nucleon, and a Finite Aluminum Slab using 
the Extended Cross Section Set 
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Figure 5.22: Scalar Neutron Flux Profile by Energy Group with a Fluorine Beam 

of 1000 MeV per Nucleon and a Finite Aluminum Slab using the 
Extended Cross Section Set 
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Figure 5.23: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 20.0 MeV (Energy Group 1) for a Semi-infinite 
Aluminum Media using the Limited Cross Section Set 



Figure 5.24: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 10.0 MeV (Energy Group 2) for a Semi-infinite 
Aluminum Media using the Limited Cross Section Set 
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Figure 5.25: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 1.0 MeV (Energy Group 3) for a Semi-infinite 
Aluminum Media using the Limited Cross Section Set 



Figure 5.26: GlT Neutron Source Profile for Energy Groups One through Four 

for a Fluorine Beam of 1000 MeV per Nucleon and a Semi-infinite 
Aluminum Media using the Limited Cross Section Set 
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Figure 5.27: GlT Neutron Source Profile for Energy Groups Five through Eight 

for a Fluorine Beam of 1000 MeV per Nucleon and a Semi-infinite 
Aluminum Media using the Limited Cross Section Set 



Figure 5.28: GlT Neutron Source Profile for Energy Groups Nine through Eleven 

or a Fluorine Beam of 1000 MeV per Nucleon and a Semi-infinite 
Aluminum Media using the Limited Cross Section Set 
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Figure 5.32: 


Angular Neutron Flux P] 
x = 200 cm for Energy G 
1000 MeV per Nucleon, z 
the Limited Cross Sectioi 
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Figure 5.33: Scalar Neutron Flux Profile by Energy Group with a Fluorine Beam 

of 1000 MeV per Nucleon and a Semi-infinite Aluminum Media 
using the Limited Cross Section Set 
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5.34: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 20.0 MeV (Energy Group 1) for a Finite Aluminum 
Slab using the Limited Cross Section Set 



Figure 5.35: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 10.0 MeV (Energy Group 2) for a Finite Aluminum 
Slab using the Limited Cross Section Set 



147 



Figure 5.36: GlT Ion Flux Profile for a Fluorine Beam of 1000 MeV per Nucleon 

at an Energy of 1.0 MeV (Energy Group 3) for a Finite Aluminum 
Slab using the Limited Cross Section Set 



Figure 5.37: GlT Neutron Source Profile for Energy Groups One through Four 

for a Fluorine Beam of 1000 MeV per Nucleon and a Finite 
Aluminum Slab using the Limited Cross Section Set 
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Figure 5.38: GlT Neutron Source Profile for Energy Groups Five through Eight 

for a Fluorine Beam of 1000 MeV per Nucleon and a Finite 
Aluminum Slab using the Limited Cross Section Set 



Figure 5.39: GlT Neutron Source Profile for Energy Groups Nine through Eleven 

for a Fluorine Beam of 1000 MeV per Nucleon and a Finite 
Aluminum Slab using the Limited Cross Section Set 
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Figure 5.40: Angular Neutron Flux Profile for x = 0 cm, x = 60 cm, and 

x = 120 cm for Energy Groups One through Three, a Fluorine 
Beam of 1000 MeV per Nucleon, and a Finite Aluminum Slab using 
the Limited Cross Section Set 
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Angular Neutron Flux Profile for x = 0 cm, x = 60 cm, and 
x — 120 cm for Energy Groups Four through Six, a Fluorine Beam 
of 1000 MeV per Nucleon, and a Finite Aluminum Slab using the 
Limited Cross Section Set 







Direction 11 


Figure 5.42: Angular Neutron Flux Profile for x = 0 cm, x = 60 cm, and 

x = 120 cm for Energy Groups Seven through Nine, a Fluorine 
Beam of 1000 MeV per Nucleon, and a Finite Aluminum Slab using 
the Limited Cross Section Set 
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Figure 5.43: Angular Neutron Flux Profile for x = 0 cm, x = 60 cm, and 

x = 120 cm for Energy Groups Ten and Eleven, a Fluorine Beam of 
1000 MeV per Nucleon, and a Finite Aluminum Slab using the 
Limited Cross Section Set 
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Figure 5.44: Scalar Neutron Flux Profile by Energy Group with a Fluorine Beam 

of 1000 MeV per Nucleon and a Finite Aluminum Slab using the 
Limited Cross Section Set 
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Figure 5.45: Magnified Neutron Source Profile of Energy Groups One through 

Four for a Fluorine Beam of 1000 MeV per Nucleon in a 
Semi-infinite Aluminum Media 



Figure 5.46: Transmitted Neutron Total Scalar Flux versus Aluminum Slab 

Thickness using the Extended Cross Section Set 
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Chapter 6 

Recommendations and Conclusions 


A multigroup, isotropic scatter, neutral particle Fn transport solver for 
semi-infinite and heterogenous slabs has been created and coupled to a closed-form 
analytical solution to the galactic cosmic ray cascade, GlT, to determine the 
behavior of low energy neutrons in the cascade. This is a first step in creating a 
neutron benchmark to be used with BRYNTRN. 

The formulation uses a closed-form representation for the neutron source in 
a basis function expansion for the neutron transport solution. Therefore, 
truncation errors are limited to the number of terms used in the expansion. These 
errors are under program control and secondary to the round-off errors caused by 
the discretization of the real numbers by the computer. The round-off error for a 
64-bit floating point representation becomes important when the incident ion 
charge is twelve or higher for the GlT algorithm and ten or higher for the Fn 
algorithm. For a 128-bit floating point representation, the charge number can be 
extended to twenty-three for the GlT algorithm. The problem manifests itself in 
the recurrence relation used to determine the partial fraction coefficients, Y^”*. 


The potential for other round-off errors can also be a problem in the Fn algorithm, 
e.g., numerous matrices are inverted, large summations are made, and numerical 
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integrations axe performed. These problems can be mitigated, but not avoided, by 
using sound and proven numerical solution techniques, such as Crout’s Algorithm 
with partial pivoting and back-substitution for solution of all matrix equations. 

The MGSLAB, MGSEMI, and FNCRIT programs created for this dissertation 
can also be used without accessing the GlT algorithm. The programs can be used 
to obtain answers to problems in themselves, or can be used as an accurate 
benchmark. Also, the Fn algorithm treats some problems that would be difficult 
for an Sn based code, in particular, deep particle penetration problems. Since the 
spatial variable is not discretized, there are no discretization errors. Therefore, the 
flux at the spatial point of interest can be determined in one calculation using the 
Fn method. In contrast, the Sn method contains inherent discretization errors 
with every spatial step in the slab until the point of interest is reached. For deep 
penetration problems with many spatial steps, the errors can accumulate until the 
resultant fluxes are meaningless. 

The major problem with the Fn algorithm is the choice of the basis 
functions used in the expansion. The closer the lower order basis functions fit the 
solution, the more accurate the Fn algorithm becomes, and the number of terms 
needed to represent the solution is reduced; thus, the potential for truncation error 
is avoided. This can be overcome in part by allowing various basis functions to be 
available to the user. The initial basis function may not solve the problem, but 
trial and error should converge on the best set of basis functions to use for the 
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particular problem being studied. 

This dissertation is the starting point for a comprehensive benchmark for 
the galactic cosmic ray cascade. The major issue that still requires addressing is 
the reformulation of the GlT closed-form solution to allow for isotopes of the ions 
being created by fragmentation; that is, the inclusion of neutrons in the 
formulation. The approximation that allows a fraction of particles with a charge of 
one to represent neutrons, the f n parameter, destroys the charge conservation of 
the multiplicities and is ad hoc. In addition, more extensive cross section sets 
should be used. The 20 MeV limitation of the Endf/B V database hinders this 
benchmark and is inconvenient for normal use. If the initial ion energy is raised 
without adding the requisite cross sections, the high energy neutrons scattered 
below 20 MeV are not considered, and invalid answers are produced as has been 
shown. 

Only the Green’s function formulation of the GlT solution is used in this 
work. This allows verification of the techniques used to solve the galactic cosmic 
ray cascade. Other formulations of the GlT solution exist for an energy distributed 
source and a composite ion beam source. The basis function expansion method 
used in the Fn algorithm can also be used on the anisotropic scattering transport 
equation. In addition, a simultaneous method or a source iteration can be 
incorporated to include flux dependent sources such as upscatter or fission. These 
utilities could be incorporated into this benchmark. 



To determine the coupled charged particle and neutron fluxes another 
method can be used. This method is based on the Sn solution to the three 
dimensional neutral particle and Spencer-Lewis charged particle transport 
equations utilizing SMART scattering (References [2] and [3]). Whether this 
method is used as a benchmark for BRYNTRN or as a replacement, it should be 
investigated when a production version of the program is released. 

The algorithm and programs presented in this dissertation generate results 
that can be verified and explained. An example and application are presented and 
the results are well understood within the physical assumptions. Therefore, this 
work, assuming the use of proper cross sections, will generate a benchmark quality 
solution to the galactic cosmic ray cascade and other transport problems of import. 
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Appendix A 

Derivation of the Analytical Solution to the Neutral Particle 

Transport Equation 

In this appendix, the steps involved in deriving the simplified Boltzmann 
equation used in the Fn method and the manipulations needed to generate the 
resultant equations used in the various programs are shown. 

The full Boltzmann transport equation is 

-- + V r + a-V v + <r(r,E) = 

v ot 

= f°dE' [ dfi'<7(r,E') f(r;ft',E' -► ft, E) $(r,n',E\i) + ( A>1 

Jo J 4tt 

+ Q(r,S7,E ,*)• 

The sections below use various atssumptions and approximations to reduce this 
equation to a form that describes the physical situation being studied and can be 
solved using analytical techniques. 

A.l Creation of a One Dimensional Transport Equation 
These steps transform the Boltzmann equation in to a one dimensional, 
multiple group, multiple region transport equation suitable for the Fn method. 

1. Assume 


Steady State. 
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• No external forces acting on the particle 

a-V v = 0. (A.2) 

• Plane geometry 

ft • r = cos 0 = p, (A. 3a) 

ft • V$(r, ft, E) = n, E), ( A.3b) 

cr(r, E)$(r, ft, E) = cr(x, E)$(x,/x, E), (A.3c) 

/°dE' [ dft' <r(r, E)f(r; ft', E' - ft, E)$(x, E') = 

Jo Ji* (A. 3d) 

HdE' / dft'<7(z,E')f(x;ft',E'^ft,E)$(x,/i,E), 

JO J\* 

Q(r,ft, E) = Q(x,/i,E). (A.3e) 

• f depends only on ft' • ft = /io 

<r(x,E') f(x;ft',E' - ft,E) = <x(x, E') f(x; E' — ► E; /i 0 ). (A.4) 

2. Expand cr(x, E') f(x; E' — * E; /x 0 ) in Legendre polynomials 

00 9 / 4-1 

<t(x, E') f(x; E' — ► E; pL 0 ) = ^ — cr/(x; E' — *• E) P/(^o), (A. 5a) 

(=o 47r 

where, 

<x;(x;E' — ► E) = 27 tJ dfiQ cr(x, E) /(x; E' — ► E; n 0 )Pi(po), (A. 5b) 

and assume isotropic scattering by truncating the series at / = 0 

j d/i 0 cr(x,E') f(x;E' — *• E;^ 0 ) = — cr 0 (x, E' -» E). (A. 6) 
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3. Expand Q(x,/i,E) in Legendre polynomials 
1=0 47r 

where, 

Q/(x,E) = 2tt J i dfi Q(x,n,E)Pi(n), 

and assume an isotropic source by truncating the series at / = 0 


Q(x,^,E) = ^Q 0 (x,E) = is(x,E). 


4. Substitute these assumptions into equation (A.l) to obtain 


d 

n— + cr(x, E) $(x,/i,E) = 

1 r°° r+i 1 

= -j Q dE' <x 0 (x; E' — ► E) J ^ dfj,' 4>(x, fi', E') + -S(x,E). 

5. Specify the standard multigroup procedure 
/ dE'(-) = £ dE'(-). 

J 0 gl = l JEgl 


Substituting this into the transport equation gives 


^Jz + ^ Z,E ^ 


${x,n,E) = 


1 ^ E 

^ f ’"^dE' <7 0 (x; E' 


1 ^ /“Li / 

= j£/ e 

Z a' = l 7E V 


E) £V*(*,A*',E') + 


+ -S(x,E). 


(A. 7a) 


(A.7b) 


(A.Sa) 


(A- 9) 


(A. 10) 


( A.l 1) 
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6. Integrate the transport equation over E G [E s , E ff _i] and define 

<t>g{x,n) = f 9 dE$(x,/x,E), 

J E, 

s‘(i) S / E '"dES(i,E), 

fE 9 _ i /Ej- i 

crg(x) / dE$(x,/x,E) = / dE <r(x, E)$(x, E). 

JEg J Eg 


The resulting set of transport equations for g = 1,2 , ,G are 

d 

“Tz + *.W 


d> 3 (x,/i) = 

1 r E 9 _i G rE-/, 

= 0 dE E / dE' cr 0 (x; E' — *• E) x 

^ JE 9 g' = 1 

X *(*,/*, E) + is'(x). 

7. Assume the scattering term can be rewritten as 

0 VE dE'$(x,^,E0 / dE<x 0 (x;E'^ E) 

2 2-i 3 < =1 [2 e 9 / 2e, 

then define 

<7g'~ g {x) /^'"‘dE'^x^E') = 

2e 9 , 

= / ''"'dE' $(x, E') / E9 "dE <r 0 (x; E - E), 

2E gl JEg 

to give 


d , . 


<f>g(x,n) — 

= \ E cr s'-s( I ) / + ^S;(x). 

Z a' = l “' -1 ^ 


(A. 12a) 
(A. 12b) 
(A. 12c) 

(A. 13) 

(A. 14a) 
(A. 14b) 


(A.15) 
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8. Assume constant cross sections across slab i 
<f>g(x,n) = 


d 

T, + < 


= I H V'g'^g / d^^ fl /(x,/) + ^S;(x). 
z 3 '=1 - 7 - 1 ^ 

9. Use these boundary conditions, where F®’ *(/z) and F^’(^) are known 
functions 


(A. 16a) 


M x i-i 1 M ) 

= KW 

/z > 0, 

(A. 16b) 

<M X «> -A*) 

= nfw 

/z > 0. 

(A. 16c) 


A. 2 Creation of a Set of Integral Equations on a Slab Boundary 
These steps create a set of complex integral equations from the transport 
equation, equations (A. 16). The complex integrals are evaluated on the real axis 
by use of the Plemelj relations. 


1. Replace /z by —fi, multiply by e , and integrate x on [ 21 , 22 ] 

Sfl f z 2 frg* 

Bg(fi,s) - a ' dx e~ > <f> g {x,-p) = 

Ji — S J z\ 

= 2s o 7 — : av( 5 ) + oT — : s 9 (s,zi,2 2 ), 


(A. 17a) 


2 fi-S 


2 fi — s 


where, 


5) — 

*q z 1 

e * 

<T Q Z 2 

) —p) — e ~ 3 <f> 3 { z 2, “/*), 

(A. 17b) 

P“gA S ) = 

/*2 

/ dxe * 1 

r 4*1 

1 d n' 

(A. 17c) 

j 


-1 


S*(s,2i,Z 2 ) 

f 2 2 

= / dx e“ 

Jzi 

■ + s;(x). 

(A.17d) 
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2. Divide by s and integrate fi on [— 1,+1] 



a 

fi-S 




3 ) 




k g <g{s) + 


+L(s)S*(s,z 1 ,z 2 ), 


where, 


^g'g{ 3 ) — Sg'g + — ^7"^L(.s), 


1 /•+» 1 

L(s) = - / d/x . 

2 7-i a — s 


3. Multiply by e 
r+i 


/ i^-^—C,( h ,s) = <r;£l£.(s,z 1 ,* 2 )A s . 11 (s) + 

■ /_1 P 5 n'=l 


"i 1 ! 


+ L(5)S*(s,z lt Z2)e * , 


where, 


C^(/i, 5) 


e » B g (n,s) 

"qi* 2-^1 ) 

~ e 1 <t>g( Z 2,-/*), 




e * 


/v( s ) 


1 r z 2 <r t g {x~z 1 ) /*+i 

- dx e > / d/x <^o'(x, — /x). 

5 7-1 


(A. 18a) 


(A. 18b) 
(A. 18c) 


(A. 19a) 


(A. 19b) 


(A. 19c) 


4. Change s = —s and /x = — /x in equation (A. 18a), and multiply by e‘ 

f +1 tx ■ G ■ 

/ <*/x D s (/x, s) = cr^ ^ J^,(s,z 1 ,z 2 )A ff - s (3) + 

^* 1 f 1 S g'zz 1 

+ L(-s)S*(-s,zi,z 2 )e' 


S' 

9*3 


(A. 20a) 
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where, 


D g {n,s) = -e * B g (-p,-s) 

g a(*3~U) 

= 9g{ z l> W ~ e ‘ M z 1»A*)» 


(A2 

- Ifdxe-^ rv 

S J z\ J — 1 

5. Define z/q to be the zero of the infinite medium dispersion relation, A 55 (s) 
(Note: if < 1 then |i/q| > 1 and real). Then L(z/q) becomes 


lk) = - 


The integral equations, equations (A. 19a) and (A. 20a), become a set of 
constraints (related to the concept of boundary conditions for differential 


equations) 


A* * / 0 


2)e 0 , 


t* “ ^0 U O a g^g 


e 0 . 


6. Use the Plemelj relation (Reference [22]) 


lim / — . 7 = p ± *>$(»? - v), 

e— '0 77 — (iy zb ze) — v 


to restrict the integrals to the real axis. If the Cauchy type integral 


xls)= r\m, 

J- 1 77 — S 


(A. 23b) 
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is evaluated using the Plemelj relations, the new formulation is 


limx^ ± ie) = X ± { 1/ ) = / d 77 1^1 ± i7r/(i/) 
o y_i rj — v 


(A. 23c) 


Apply these rules of integration to the integral equations to obtain 


y+i 7 ] 

f d rj C g (rj,u) ± invC g (v,v) = 

J - 1 77 — 1/ 

G 


= ^ + e“^ i S;(i/,z 1 ,2: 2 )L ± (£/), 

g' = l 


(A. 24 a) 


/* + l fj 

f drj Dg(r/,u) ± iirvD g (v,v) = 

- 7 - 1 v ~Q a ,^ (A. 24b) 

= 0-3 H J ii'( I/ > z i» z a) + e“^ r "S*(-t/,r 1 ,2 2 )L ± (i/), 

3 ' = 1 


where, 


* f. 


(A. 24 c) 


Aj f (*) = A^) ± “ 


2CT 3 


(A. 24 d) 


z^cr 1 , r+i 1 


2 ( 7 * 3-1 T] — 1/ 


2 <t> 


1 - 1/ 


1 + 7 / 


,(A. 24 e) 


(note: all integration variables have been changed to 77 and u £ [0, 1] U t/q) 


7 . Eliminate 21, z 2 ) and zi,z 2 ) by adding and subtracting the 


positive and negative branches of the equations and combining, there results 
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+ 


v 


C 9 {ri y u) 


ale * 

— ^7 — S*(v, 

ua % 

9~*9 

a g—g y' i 

— 2- <V-s 

3 3=1 

s'#s 


a g-g 

21 , 22 ) + 

^("1 2l,2 2 ), 


*0 = 


(A. 25a) 


/ <ty — —Dg{v,v)--r-KM)D 9 {'', 1 ') = 

7-1 Tj — u ^ 


9^9 


cr o z 2 


<^ € * 


3 

G 


S*(-t/, 21 , 22 ) + 


(A. 25b) 


+ -rr E j;;«(^,2i,2 2 ). 


* <?'=i 

9^9 


8. Rewrite the integral equations in terms of the fluxes by changing the 
integration variables so they are evaluated on the interval [0, 1] to obtain 

yi 77 yi r] 

T di7 — <£ 3 (x,-_i,-?7) + d V — — <^ 3 (^- 1 ^)- 

Jo 7 ] — 1/ Jo 7 } + V 

&*g \ 7] 7] 

-e - d ?7 <f> g (xi,-ii)+ dTi——<j) g (xi,r}) - 

[SO TJ — 1/ JO T] + 1 / 


<7 9^9 


X 'ggi y ) 


A’„ 


<j> g (xi- u - v) - <t>g(xi, -u)e » 


(A. 26a) 


= —ir~ £ <7 s'-*s I ss , ( i/ , x *-i> x ») - 

9^9 g’= 1 

g’^g 

*g*i - 1 




^s-s 3 


1 1 ) 1 
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-fdT)—2—<t> g (Xi,T))+ f dr] —^—<t>g(xi, -T}) - 
Jo 7) — V Jo 7] + V 

£ dr) —?—<t> g {xi- U Ti) +£ di] —^—</> g (xi. u -t]) 


*? 

■t » 


<t>g(Xx,v) ~ 4>g{ x i-\i u ) e u 


2a' 

- 

a 9~9 . 

= --T" £ ag'^A^Xi-UXi) - 

3^9 g'=\ 

9^9 


A* 


(A. 26b) 


frje * 


ua 




9-+g 


where, z\ and z 2 are defined as the slab boundaries, x,_i and x,- and 
= < 7 * (x t — Xi_i) is the dimensionless slab width. 


These integral equations are in the form of an inhomogeneous Fredholm equation. 
Various methods of solution are available. The one chosen for this work is a basis 
function expansion and a collocation method for determination of the expansion 
coefficients. 


A. 3 Application of the Fn Approximation to the Slab Boundary 

Integral Equations 

The method used to solve equations (A. 26) is a basis function expansion 
called the Fn approximation first developed by C. E. Siewert. The basis functions 
can be any set of functions; however, orthogonal functions on the interval [0, 1] 
create matrices that are not ill-conditioned as with the matrices generated by 
non-orthogonal functions. 

Substitute these boundary conditions and Fn approximations into the 



169 


integral equations 


4>g( x i-U ") = F£‘(t'), 


(A. 27a) 

0 fl (x, -,-!/) = F §*(»/), 


(A. 27b) 

<t>g(xi-i ,v) = F^‘(i/)e - "^ + 

<j' JV-1 

2~’ E w-m, 

^ a=0 

(A. 28a) 

<f> g {xi,-v) = Ff l ' , (i/)e""^ + 

<7* A^— 1 

-£r E <•'«-). 

(A. 28b) 


g a=0 


to generate a related set of singular integral equations that are solved for the 
expansion coefficients using a collocation procedure. The resultant integral 


equations are 


iV-l 

£ 

a=0 L 


_±i 




X «-l) + 


(A. 29a) 


+ — — + — — Sl^X.-^X,), 


9—9 


9—9 


N - 1 


ot — 0 L 


+ a£' i A»(i/)e* 


= R2* (*/,*,) + 


1 


1 


(A. 29b) 


+ — — T2j(j/,x,_ 1 ,x I ) + — — S2^(t', x t _i, x t ), 


9^9 


9-^9 


where, 


KM = 

~ '’ITT ~f ^7 

2a' g Jo T) — i/ 

= J‘ivv{ Fr(,)C(A;,^,,) + Ff (,)S(A>,,)] , 


(A. 29c) 

(A.29d) 

(A.29e) 
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R2‘ (i/,x,) = J\v[ F£' , '(»;)C(A; i i/,iy) + F&‘'(j?)S(A;,i/,7?)] , 

(A.29f) 

e * — e * 

• 

(A.29g) 

1 _ e - a i!i+;] 
S( A i. „ ' 

(A.29h) 

G 

Tl’Cl/.Xi-!,*,-) = (T\ Y, 

9 \ 

(A.29i) 

9 =1 

g'*g 


G 

T2^(i/, X,_l, Xj) = (Tg ^ (Tgt^gjggf^^Xi- l,X t ), 

(A.29J) 

3 /= l 

9^9 


SI* (|/,X, ■_!.*<) = — / dze * S g (z), 

V Jx % - 1 

(A. 29k) 

(7*- /* r « *o(*»“*> , . 

S2;(i/,Xi_!,x.) = — / dze~ - Sj(z). 

V Jx % - 1 

(A. 291) 


These equations are changed into matrix equations using the collocation procedure 
and the expansion coefficients a^’ 1 and b^’ 1 are determined using a matrix 
decomposition scheme as described in Chapter 3. 

A. 3.1 Singularities Encountered in the Fn Equations 

The integral term the exponential term C(A^, v,tj), and the source 

terms Sl^(t/, and have singularities which must be carefully 

handled. 

For the B|(i/) term, the original definition is 



B s » = \ aa {v)M») ~ 


(A. 30a) 
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where, 


cr x 1 / /■+ 1 1 

*«(■') = 1 + ~§h £, *' 


(A. 30b) 


To regularize this term, if 

/ dT ? r~ 

y-i 77 — 1 / 

is added to the second term and subtracted from the first term, then B£(i/) 


becomes 


B*M = A- - b;»(i/), 


(A. 31a) 


where, 


A - M = ! _ $=£ j\ J- = l + S=» 

' 2 <r‘ Jo tj + 1 / l<rj 

b-»(„) = <=» /a, gtahl ~ 

Q v ’ 2a' Jo n-v 


— In — Z 


l + «/P 


(A. 31b) 


(A. 31c) 


When v = 77 , L’Hospital’s rule is used to find a new integrand for B‘ 9 (i/), or 


<t« ,1 v -T-i}) a {v) + tM") n = v 


(A.31d) 


T) — v 


T)^U. 


For the C(A*i/, 77 ) term, when u = 77 , use L’Hospital’s rule to obtain 


A* a* 

C(A‘ 1 ^) - 


(A. 32) 


The source terms encounter a singularity as v — ► 0. This problem is solved 
by structuring the terms to conform to a delta function definition 


[*i 

Slg(j/,x,_i,x,) = a' g / dz ■ 


_ 1 ) 


sit*), 
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(A. 33b) 


As u — ► 0, the functions in the brackets behave as the delta functions 
S(cr' g (z — x,_i)) and <5(<r‘(x,- — z)). When the delta functions are placed in the 
source terms and evaluated, there results 

r 

S’(xi_i) u = 0 

Slg( I/ ) Xj-ijX,-) = { v _ (A 


s* (x,-_i) 

v = 0 

cr\ f x i 

-*/ dz e 2 » S^(z 

V Jxi- X 

I "7^0, 


i/ = 0 

<ri r x * ■ 

— I dze » Sg(z) 

V Jxi-i 

v± 0. 


S2‘ (i/, x,_i,x,) = 


A. 3. 2 Post Processor 

With the regularized B£(y) term, the Fn approximations can be rewritten 
to achieve a faster convergence rate (fewer number of expansion terms required to 
converge). These equations comprise the post processor and are used to determine 
the final values of the flux. 

If the non-singular expression for B %({*) is substituted into the integral 
equations, equations (A. 29), the results are 

lf*;'(W.W - BJM) + = 

“=° L J (A. 34a) 

= Ri;(/i,Xf_i) + — Tl^/^X,-.!,!,) + — Slj(iU,Xj_i,Xj), 


t, b ’.'' 


= R2j(^,Xi) + -i— T2’(/i,x,_ 1 ,x,) + — — S2‘(/i,x i _ 1 ,x,-). 
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These equations can be solve for the basis functions to obtain 

E «4 'iMaO = TTTvf ( E + 

a=0 l La=0 J 

+Ri;(/i,X t - 1 ) + — — Tl* (/l, *,•-!, «i) + — f — Sl^(/X, Xi-1, X,) 1 , 

^—<7 a 9->9 ) 


iV-1 


E bi'V.(/i) = t: 


or=0 


a »(a> 


E b i’‘B; s W - 


a=0 


+ 


+R2*(/i,x,) + — — T2'(/i,x,_i,x 1 ) + E— S2’(/i,x i _i ) Xi)l . 

a ^a— a ) 

These equations are substituted into the Fn approximations, 

equations (A. 28), to obtain the post processor 

M x i- u -/*) = E ^ { <7 J-» R liO l >x.-- 1 )+ 


iV-1 


+<T a-a E 


a=0 


+ 


a^*B^( M ) - b^A^( /i )e- * 

+Slj(/i, X,_ x , X.) + Tl^(/i, Xi_i, X,)} , 

= FfOOe"^ + 2 ^ ' a* ( m ) { cr fl~ t, g R2 a(^’ I, )+ 


iv-i r 


+<7 a-a E 


a-0 


.ek 


K''KW - 


+ 


+S2’0z,Xj_i,x,-) + T2^(/i, Xj_i,x,)} . 


A. 3. 3 Scattering Terms 

In this section, the expressions for x,_i, x,) and x,_x 

equations 


Ti;(i/,x,-_i,x,) = £ v^A^Xi-uXi), 

g '= i 

9*±9 

G 

T2^(z/, Xj_ x , x t ) = <7g E 

a'=i 

9^9 


(A. 35a) 


(A. 35b) 


(A. 36a) 


(A. 36b) 


x,) in the 


(A. 37a) 


(A. 37b) 
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for the scattering terms axe derived. 

Since the particle group transfer functions x,_i, x,) and 

x,_i, x,) axe dependent on the spatial integral of the scalar flux for all other 
groups, the scalar flux must be known at all points in the slab. The numerical 
scheme used to solve the Fn equations dictates that the boundary values must be 
determined before the interior values. These steps show how the transfer functions 
axe reformulated in terms of the other group angular fluxes at the boundary. 




3. Multiply equation (A. 39) by e ** and equation (A. 40) by e ** to obtain 


4>A Z 2>A«) - 4>g'( z ufi) e ~ " 


n 

G 


s *A-H, z U z 2) + 


1 f z 2 

+ - z2 a g"~g' / dl 4g"( x ) 
P 3 "= l Jz ' 


e m 


(A. 41a) 




e * 
G 


S*/(/i,Zl,22) + 


1 Jf, y-22 

+ - <»-«,' / dx 0 3 »(x) 

t* g" = l Jzx 


e ** 


(A. 41b) 


4. Transform the exponential in the scalar flux integral to be a function of 5 

^ i 

and not by letting /i = -^r£ = s gg'£ < 1 so that £ £ [0, to obtain 


47 ( z 2 “ Z 1 ) 
a „/£ 




_^2 
e 3 99 /€ ^ 

SgA ^ 


S^(-5 w ^,Z 1? Z 2 ) + 


*’,(*2“*) 


+ Y1 A-’-g 1 / d *<M*) e >9,,< > 

5 55'S gtt-l Jz \ 


(A. 42a) 


<7 g/( z 2“ z l ) 
” * r,„> £ 


It $gg'£) $g , { z 2y $gg J £)^ 99 


n f Z l 


e 


Sgg'£ 


S m M 9 gi^Z U Z 2 ) + 


(A. 42b) 
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5. Remember — = <x* so that the transfer functions are rewritten as 


1 fZl 2 - j ) 

- dxe « <f>g"{x) = J 


1 r«2 

- I dxe < V( x ) = ^"(f^i^a)- 
s •'*1 


(A. 43a) 
(A. 43b) 


Substitution into the transport equation solution gives 


4>9'{z2,s g9 '0 - <f>g‘(zus gg .()e a 99‘ 


2 ”*l) 


<* q , z T. 

e a 99 fi 

S 3 A la 

~i G 


S 3 ' ( s gg'£) z i, Z2) + 




‘V 3''=! 


(A. 44a) 


< 7 g/( x 2 -‘ E l) 


5 55 ; 0 S 99 , ^) € ' ,j9 


n t Z l 


e 

s gg f £ 

G 


S*A S 99'tZl,Z 2 ) + 


(A. 44b) 


+ -f E 

<V 3"=1 

6. Again, let zi and z 2 be the slab boundaries x;_i and x,*, then the above 
equations reduce to 


*nn't 


'{Xi,Sgg '0 <f>g*{Xi—\ , Sggt^)t 9? 


J ,f 

e « 


-f ^i-li 2?|) + 

5 35's 


*99 

r x G 


+ ~T H < 7 3 "-g' J 33 "(^ X «-X’ X >) ? 


‘V a"=i 


(A. 45a) 
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*L 

( t*g t {Xi—h~3gg t £) — <j>gt{X{, ~Sgg'£)e 99 


y x »- 1 

e 


(A. 45b) 


> ^t-l 5 ^t) + 

^ss's 

+ £ E 

*V g"—l 

7. Apply the boundary conditions and the Fn approximation for group </, and 
solve for !*’,(£, x,_i, x t ) and 




N-\ 


' 9 a=0 


. Sl 5 /(s 55 '£ , X{- 1 , x<) 

r 5 

G 

E X t-1? ^t)j 


(A. 46a) 


*?"=! 


^ Ar " 1 




c*=0 


,• S2 5 /(s 55 /^, £*-ii ^i) 

G 

a ’g"-'g‘J gg"{£l X i-ll X i)- 


(A. 46b) 


9 " = 1 

p'Vs' 


These equations can be used directly in Tl^(/i, z^-i, x,-) and T2‘ (/*, x,-_i, x,-) 
for <f £ [0, jiy] as long as is not equal to zero. 

8. For if ^ [0, — ], i> > 1, the results from equations (A. 19a) and (A. 20a) for 

s 99' 

group g' can be used because they no longer generate singular integrals 

f+1 n 

/ <ty— — C a .(i/,j/) = <r‘, X, Ig! s »( 1/ ^ 1 ’ 22 )A9"3’( ty ) + 

• / -i V v » =1 


+ L(i/)S*,(i/,Zi, z 2 )e‘ 


v 1 ' 


(A. 47a) 
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f+l n . J® . 

/ drj Dgi(t), v) = ^ XMsV(^ ^1 ’ 22 ) A s'v( I/ ) + 

•'-l »7 - 9 ^ 


S' 

/*2 


+ L(i/)S‘,(-j/,z 1 ,z 2 )c - , 


where, 


C 5 /(r/,i/) = - e - <t>g’{Z 2 , —Tf), 


gg(»2~*l) 


Dg'iv,^) = 4>g'{z2,Tl) - e - <j> g .(z U T]), 


1 [ Z7 «' f 

I*; s »(i/, Zi,z 2 ) = - j die' - y_ t dr/ <£ 3 »(x, -t?), 

1 z- 22 g a<(*a- x > y+l 

J*l„(i/,Zi, z 2 ) = - dxe * / d77 (f> g »(x, -rj), 

* V J Z\ v/ — 1 

a s v( j/ ) = Vv + ^ s r 5 l(v), 

°V 

1 /+ 1 1 

L(") = 9 / » 

S*,(i/,z 1 ,z 2 ) = f dx e~^rS' g ,(x). 

J Z\ 


(A. 47b) 

(A. 47c) 
(A.47d) 
(A.47e) 
( A.47f) 
(A.47g) 
(A.47h) 
(A.47i) 


9. Rewrite the integral equations in terms of the flux values by splitting the 
integrals over the interval 77 £ [ 0 , 1 ] to obtain 

[l „ yl „ 

/ <M 2 1 >~ 7 ?) - e - / dT? <j>g'{z 2 ,--q) 

Jo T) — V Jo T) — V 


+ 


+ 


Jo*” 


T) + V 




— e 


g ^( z 2~ *1 ) r 1 


/ d»? — <t>g>(z 2 ,T}) = 
Jo T) + V 




sr"=i 


+ L(i/)S‘<(j/, zi,z 2 )e - ^ - , 


(A. 48a) 
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fl T) rl n 

/ <1 t) <t> S '(z2,v) - e v / dy + 

JO Tf — U Jo T] — V 

+ / d»? —7— </>A^2,-y) ~ 

Jo T] + V 

-e * /dr? — - — = 

JO T] + 1/ 


— e * 


= ^3"s'( t/ )JsV'( i/ ’ 2l > Z2 ) + 


+ L(t/)S*,(-t/,Zi,z 2 )e' 


10. Again, let Zj = x,_i and z 2 = £;, apply the boundary conditions and the Fm 
approximation, and substitute in known functions 


£ - bJ'AJ (i*)e " + Ri;, (*,*,-.) - 

a=0 - 

- L(i/)S’,(i/,x,_i,x,)e * = 

G 


53 + R2;,(i/,x,) - 


- L(i/)S*,(-i/, x,_i,x,)e = 

G 

= 53 A s - v (i/)j*: s „(i/,x,_i,x t 


11. Let // = = Sggt £ > 1, then 


<7 ‘,I!L„(l/,X t -i,X,) = <!!'»(£, x.-i,*;), 




(A. 50b) 
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12. Solve for ^gg > (£} *e»— i f ^*) 3 >nd. i? *^*) 


X i-1 > ^i) “ x i-l ) 

- *<-!,*■) + 


AT-1 

+ £ 

cr=0 


A 1 , 


a> v A ;'(-<»'{) " e“-..' , b;'- i Aj'( S „.« 


1> •*■»)> 

s"=l 

3"#9' 


(A. 51a) 


^g^-9'9 1 (£ ) ®t— 1 > ®« ) ^-2g/(s a g'£, X,) 

- :j, + 


AT— 1 

+ £ 

a=0 






G 

~ &g y! ^S"3'( 5 S3'£)Jg3"(£> ®i— 1> ^i)’ 
s"=l 

9 "* 9 ' 


(A. 51b) 


These equations can be used directly in Tl^(/x, x,-_i, x,) and T2 g (fi, x,_i, x,) 
for £ ^ [ 0 , — ] as long as A s / s /(s ss /£) is not equal to zero. 

3 99' 


A. 3. 4 Problems Encountered in the Scattering Terms 

As noted in Section A. 3. 3, the terms for x;_j, x.) and x,_i, x,) 

have singularities at 


• For £ £ [0, with cr' g ,_ g , = 0. 

• For £ £ [0, 7 ^ 7 ] with A g >g>{s gg >Z) = 0 or s gg ^ = t/£. 
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The steps below determine the expressions that axe used in Tlj(/z, 
and T2*(/i,x,-_i,Xi) for the points indicated above. 

1. For £ 6 [0, jhj-] with = 0, staxt with the singular equations 

= ° J ^T £ a-a’^or (Sgg'O - 


9 a=0 


: Si gi { s gg'£i 3 Ci — 1 7 3 f «) 

a 9 

G 

— (T 9"-~9‘ ^33" (^’ X i-li I ')) 


s"=i 

a"*g' 


N - 1 




7 5 £>=0 


• S2„/(Sjg'£, X,_l , x,) 

a 9 

G 

- a ' g "^3‘^ gg "{^ X ^- UX t ). 


g"=\ 

9"*9‘ 


Since cr‘, , = 0, the equations above reduce to 


) ' (X g"—g'lgg"{£l X i—ll X i) — { ^ lg'^jg'^, X,_ J , X, ) , 


P" = 1 

9"*9‘ 


&g"—.gl ^gg"(^ 1 — 1 > • J '») ~ , S2 9 /(Sgg<^, X,_l, X,). 


3 " = 1 


2. Since <7^ s , — 0, these terms become 


Ag/g'(Sggl£) 1, 


^■9"9'( 3 99'£) ~~ ( S 99'^)^-‘( S 99 I 0 i 9 > 




®i— 1 j ■£«) 


(A. 52a) 


(A. 52b) 


(A. 53a) 


(A. 53b) 


(A. 54a) 


(A. 54b) 



182 



These equations are used directly in T1‘ (/i, x,_i, x,-) and T2' g ( fi, x,_i, x,) for 
£ € [0, —1 when < 7 *. , is zero. 

^ L S gg* 1 9 9 
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4. For £ £ [0, j^-] with A g ' g '{s gg '() = 0, or s gg >£ = i/g, start with the singular 
equations 


“ A lV h„.() 

- ^yvt isi •, . m , + 

r a 9' 

N-\ ^ 

+ E 4' a £ (-<«■() - - 

or=0 




:r «) "**"* a / X{) 

i'a'a'ySgg'S) 


- M)l f a "' <) st,M,x l . l ,x l ) + 


+ E - e _ -..-'aJ , ’ i AS'( Sss .{) ! 


0*0 (£ 5 ^t-l 7 ^t) f * 


5. Take the limit as s gg >t —* through the use of L’Hospital’s Rule 


lim = 1 7 — 7 7TT7~ 71 1 

■ , 39 , ^ _ " / 0 d a £ Ap'g'fSjj/^) ds 3fl <^ 

) Cl , t _ t _ _ \ , 

■ ol 5 i(S ss '^, Xi_x, X,) + 

<v 

,'V-l A '„< 

+ E a a’*Aa (~ s a f '0 - e'^b^A'WO 


{RlyCs^^Xi-!) - 


Xi-1 , X t ) > , 
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lim 


°g^7A^ Xi - 1,X< ) ~ d ^ / ,£)d 5 ,£ { R2 3'( 5 33'^ 

di ,(W\Vsj a, 5 P 5 c 






N-\ 

+ E 

a=0 






G g XI X,_J , Xj) 


s"=l 
a" *9’ 


> . 


6. These differentials are evaluated term by term 


d-Sgg'£ 


-AiM = %?j( > urhr^ M- 


A* , 


dSaa # £ 


e A£'(<wO = e '*** 


(V S gg '£) 2 

d 

d$gg*£ 


*IM + 


+ (»»■<)' AS(3 "' 0 


d . , * \ Gg"~tg' 

T" 7^9"9'( S 99'0 — n i 
dSgg'Q 2 &g> 


In 


f s gg'£ 


1 d" s gg'£ 


+ 


2 (fgA) 

(^) 2 -iJ 


d5 ^Rl^(s^<f,Xi_i) = jf drj ?7 [FS[**(i|) X 


x :- -- . C(Ah,T/,3 3J .0 + F^i?) d g S(Aj,,iy,s w ^) 




d^ 55 ^ 


di — ^R2;,(s 33 .{,ij) = [f{''(>() x 


x jd-c(A;„,,3 3S .o + 

dSj fl <^ ds 3P <f 


K> 


ddpg'^ 


C(A g,,ri,Sgg>Z) = 


C(Ag',T 1 ,Sgg'£) (j ff ) » C 


T) ~ Sgg'Z 


! ®» ) 

(A. 58b) 

(A. 59a) 
(A. 59b) 
(A. 59c) 
(A.59d) 

(A.59e) 


(A.59f) 
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dSggt^ 


S(A*„7M^0 = i-T [S(A;„ 77,^0 + 


+ 




(^■O 1 


d ✓ x,*_ i, X|) 


d-Soo'^ 


( ,S 55 / 0L( 5 55 / 0' 




( ( d SI , Xi-i, x,*) 

= (*,0.OL(»„-O 3^ • 

+ Lfa,.{) S1 *' (, "'^’ Ii ''’— 


a ' 3 ‘ 


+ 


a * 


+ 


+ 


s gg'€ Slgt{Sggt£, Xt-1 J Xj) 

MY - 1 4 


d , ^ ^ S2y X,_i , X, ) 


dSgg/£ 


(Sgg'O^iSgg' 0“ 




_ / >. u / c v d S2 g j(<s Jg /£, x »-i> x t) 

— ( 5 gg's) Mbps's) 


d“ 


d<s^/f 

^J-Ii ^i) 


V 


+ 


S2^/(5 5Sr /^, x,_i, x,) 

(MY - 1 


+ 


ds gg '£ 


h.g"g'(Sgg'£Ygg"(Sgg t £') x i-h x i) — 


~ ^gg ,, ( S 99 , £i x *-l i X i ) , c ^g n 9 , ( S 99' 0 d" 




+ hgltg^Sggf^) lg gi ft ( Sgg ' ^ , X { -_ 1 , X { ) . 

d^gg's 


d<Sgg'£ 


^g ,, g , ( s gg , V)d ggtti^Sggi^ X t _i , X,) — 


— dgg" ( s gg'£i x i-l 7 x t) j ^ d" 


dSgg*£ 


d” j „ t ’ ^t-l ? ^i)- 


_d_ 

d^gg 1 £ 


(A.59g) 


(A.59h) 


( A.59i) 


(A.59j) 


(A. 59k) 
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d Sl„/(.s s3 '£, Xj_i, i|) 






1 f*i 


1) 


I dzSj,(z)e '»»'* x 

Ji,_i 


(5 33 '0 3 

[ * ( <7 s'( z -*•-») -Sga’t) 

d Xi_i , 3?i) 


* 9 A = 0 




d Sgg'£> 


*9' 


If d.Sj.We-^x 

j'£) J Vr,_i 


( 5 55 # 0 

( X (<^( X > “ Z ) “ 5 5S^) 


$gg f £ 0 


Sgg't ^ °- 


(A. 591) 


(A. 59m) 


The terms Vg t g ,, i s gg , £'> x *-i i £«) an d ^i-ij^i) can n °l be 

evaluated except if there is a limit of strict down scatter and g f only being 1. Then 
the terms do not need to be determined. This calculation is used only to verify the 
output in Reference [12] and is only contained in the original non-production 
version of Mgslab. 


A. 4 Creation of a Set of Integral Equations for Interior Slab Points 
These steps determine the matrix equations and post processor for the 
interior flux. Rewrite the integral equations, equations (A.25), in terms of the 
fluxes by changing the integration variables so that they are evaluated on the 
interval [0, 1]. Define z x and z 2 as x and £; for equation (A. 25a) and as x t - X and x 
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for equation (A. 25b) to obtain 


/ d»7 — — <f> g (x, -Tj) + [ d Tj - 
Jo Tj — V Jo 77 + v 


V 

— e 


*«(*.-*) 


5—5 


1 J , + 


/> 


7] — V JO 7} + V 

4>g{x,-v) - <t>g(Xi,-v)e 


a 9~9 g '=l 
9*^9 


E 4- s I i'( i/ ’ x ’ x -) - 


a g e " 


V<J 


s m g (v,x,x,), 


5—5 


/ d?7 — — <j>g{x,ri) + / dr/ — <j>g(x,-Ti) - 

Jo T) — v Jo r] + 1/ 

<^(*-*<-1 f t 1 

- e_ ’ [to' 


V + 

nM x i- hi) , ( 1 a 






5—5 


7} — V 

<t>g{* ,u) - <t>g{Xi-i,l/)e 


TJ + 19 




<7-5 <,'=1 
5^5 

<tLj: 




i/<7 




5—5 


(A. 60a) 


(A. 60b) 


A. 5 Application of the Fn Approximation to the Interior Slab Integral Equations 
The same method used on the boundary integral equations is used on the 
interior equations using different expansion coefficients. The boundary conditions 
and Fn approximations are 


<f>g{Xi- 1,V) = Ffl'(v), 


(A. 61a) 


<f> g (x t ,-v) = F&‘("), 


(A. 61b) 
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i JV-1 


<i>g{x,-v) = Fr'(i/)c - + -2=2. J2 C™1p a {v), 


(A. 62a) 


' 9 a=0 


M*>») = + ^2 0 ^). 


'y a=0 


Upon substitution of the boundaxy conditions and the Fn approximations into the 
integral equations, a set of equations is produced that can be solved for the 
expansion coefficients 

- d"A»»] = Rl' (►,*) + 

Q=0 (A. 63a) 

+ -r— TV g {i/,x,Xi) + -j— Slg(i/,i,x,), 

*9-9 *9-9 

£[d"BiM - cJ j a;(i/)] = R2;(j-,x) + 

«=° (A. 63b) 

+ — — T2'(i/,x,_i,x) + — — S2j(i/,x i _ 1 ,x), 

/T* /T* J 


where, 


RljK*) = [FR , ( , ?)C(cr‘(Xi — x), V, T}) + 

+ e" S< ‘ F['‘(t/)S(^(x, - x),i/,r/) , 

R2*(i/,x) = J o ^VV [Fl‘(^)C(<t^(x - x,_i),i/,t/) + 
+ e"^ F^‘(t 7 )S(<t'(x - x<_ 1 ),i/,77)| , 


T1 \(v,x,Xi) = a\ ^2 <7^_ a I”,(j/,x,x t ), 


T2 s (j/, x,_l,x) — <T 5 Ggt^g X i-l i X )i 


(A.63f) 
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o' r*i 

l'(v,x,Xi) = — / dze 
V J X 

l' g {v, x,-_i,x) = — / d 


slW, 


(7* r 37 «q(*»-«) 

— I dze - S‘(z), 


(A.63g) 


(A.63h) 


These equations are changed into matrix equations using the collocation procedure 
and the expansion coefficients c 9 j and d£ ,J are determined using a matrix 
decomposition scheme as described in Chapter 3. 


A. 5.1 Post Processor 

To accelerate the convergence rate of the Fn approximations, 
equations(A.62), the regularized B£(/x) from equation (A. 31) is used in the interior 
integral equations, equations (A. 63), to obtain 


KM) £ - £ [c^Byoo + d^A »(/!)) = 

a =0 a — 0 

= R1 gift, x) + -4 — Tl^(/i, x, X.) + -J— -SI’ (/I, X, X,), 


(A. 64a) 


KM) £ WMn) - £ [d 9 a : &M + c"Ai(,o] = 

ar=0 Of — 0 

= R2‘(/i,x) + -4 — T2 ^(/x, Xi_i , x) + -4 —S2' g (/i,xi-i,x) 


(A. 64b) 


These equations are solve for J2a=o c a^a(M) an( i I2a=o d®’-*0 o (//) and 
substituted into the Fn approximations, equations (A. 62) to obtain 


4 > 3 { x ,- h ) = F&‘(/i)e {a;_ g Ri;(/i,x) + 

C S J B?W+ ( - bfe- — A'M + 

a=0 L V / 

+Sl’(/i,x,x,) + Tl^(/i, x, x t )| , 


" A*(/i) + 


(A. 65a) 
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M x *l*) = F h(l*) e ^ ^ + 2cr » A* '(/ O { <7 ^ gR2 ^ AI,g) + 


AT-1 


E 


of=0 L 


dS i B7(/0+ (c^-a^e 


^(x-x^.j) 


A£(p) 


+ 


(A. 65b) 


+S2^(/i, i,_i, x) + T2g(fi, Xi-i, x)} . 


A. 5. 2 Scattering Terms 

Using the same methodology as for the boundary, derive expressions for 


9^9 

(A. 66a) 

G 

a l E 

(A. 66b) 


<7'=1 

9^9 


1. From equations (A. 44), let z\ and z 2 equal x t _i and x for equation (A. 44a) 
and x and x t for equation (A. 44b) to obtain 




<f>gt(x , Sggt£) (f> g ' ( X { «. 1 , Sgg> <f ) £ 99 


e *99'* 


c Sgg*£ ? ^t-1 ? 3?) 4“ 

Sgg’Z 


*99 

r* C? 






(A. 67a) 




^/(x, 'Sc/j'O s gg'£) e 99 


e*99'* 


c X, > ®*) "t" 

5 55 / S 


*99 




fl"=i 


(A. 67b) 
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2. Apply the boundary conditions and the Fn approximations for group g and 


solve for !*’,(£, x,x,) and J*‘,(£,x,_!,x) 


-i N - 1 




'5 Cr=0 

G 


,■ Slg'(Sgg'£, X, X t ) ^ X ’ X *)> 


(A. 68a) 


s"=i 


i jv-i 

<■ 7 o'-in' 


= -£ ~7 2 " E d(’ J ip a (s gg .£) - 

^3 o=0 

1 a 
~ jS2gi(s gg i£, x t _x , x) — o’ g n_ g i x,_i , x). 


(A. 68b) 


s"=i 

j'W 


These equations are used directly in Tl^(/z,x,x,) and T2^(/z, x,_i, x) for 
£ E [0, jL_] as long as <? x g ,^ g , is not equal to zero. 

3. For £ [0, — ], start with equations (A. 48) for group g' . Let z\ and z-i equal 

3 99 ' 

x and x, for equation (A. 48a) and x t _i and x for equation (A. 48b). Apply 
the boundary conditions and the Fn approximations 


AT-l 

E 

a= 0 


<4' J AJ'(-I/) + (d£ J - I a;» 


+ 


Rlj,(t/,x) - L(i/)S m gl {v,x,Xi)e^ 

G 

= a \- E A «v(*') I Jv'( l '» x »*i)» 

3"=1 


(A. 69a) 


N-l 

E 

Qt=0 


d^A*'(-i/) + fc^-a^e- V< ^ X ‘" ,) ) Ai(u) 


+ 


+ R2’,(i/,x) - L(j/)S*/(-t/,x,_i,x)e‘ 

a 

= CTg' E A 3"3'( i/ ) J 3V , ( l/ ’ X ‘'- 1 ’ X )- 


3 " = 1 


(A. 69b) 
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<7* t 

4. Let v = = s gg <£ > 1, then 

£r’J*; a »(i/,x,x,) = cr‘I* a „(^,x,x,), 

&g> Jg'gf'fal ®i— 1 ) ®) = (£> **'•’- 1 » ® ) * 

5. Solve for I* a ,(£,x,Xi) and J* a /(£, Xj_i, x) to obtain 


( 7gh-g'g'{s g g'£)1- g gi(.£i x, Xj ) — Rl 3 /(5 a3 <^, x) 


iV -1 




+ E [< J A;'(-s„.fl + 
+ 


Qf =0 

( 


* b <* i -*> 
/ * . 


d^._b^*e I (Sgg'O 

G 


^ i ^ 5 // 5 , ( , 5 55 /( ^ • r *)? 


<?"=! 


(J g^9 t 9 / {Sgg'£) ^ gg f (£1 *^»-l ? **0 — (s 55 '£ ? x) 

( 5 W # 0 ^( 5 W # 0 Qrti / A \ . 

v 99*^1 ^ t — 1 ? i 


iV -1 


<’;■ 


+ E <£■**£(-*«■() + 


c *=0 


/ *— 


+ c ^- a^e A '(*„«0 


X ,_ x , x). 

sr"=i 


These equations are used directly in Tl^(/x,x ? x;) and T2^(/x, x,_i, x) for 
£ & [0, — ^ — ] as long as A $v(«w£) is not equal to zero. 

■* aa , 


(A. 70a) 
(A. 70b) 


(A. 71a) 


(A. 71b) 
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A. 6 Boundary Conditions for Beam and Isotropic Sources 
In this section, the boundary conditions are determined for multiple slabs 
for a beam or an isotropic source incident on the leftmost face. The slabs are 
connected through the Fl‘(^) and Fr‘(/j) terms in Rl‘(/i,x) and R2^(/x,i) and 
the post processor. 

To obtain the proper relations for a beam source, define the boundary 
conditions at the slab boundaries as: 

Leftmost slab - Beam Source on left face 


Ffr'M = (A. 72a) 

Fft'oo = *,(*1.-/0- ( A ' 72b ) 

Rightmost slab - Vacuum on right face 

f l’ NS (^) = ^s( x ns-i,a0 + Sq6~^ ^*=i A ^(/i - /ig), (A. 73a) 

Fr NS (/^) = 0. (A. 73b) 

Non-boundary slabs - Collided plus uncollided flux on both faces 

F £’*(/*) = <A a (x,_ x ,/x) + - n 3 0 ), (A. 74a) 


Fr ‘(/0 = 


(A. 74b) 
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A pattern is set up when F^'(p) and Fr‘(^) are substituted into Rl^(/z,z) 
and R2j(/i,x) which is 

RV g {p,x) = ldriri[C(<T , g (xi-x) 1 Ti,n)<j>g{xi,—r}) + 

J 0 


gq(x~J f -.l) 

+ e - S{a' {xi - x),r},n)(f> g (xi- U T]) 


+ 


+ a si swifr-x),^). 


R2' t (n,x) = J dr) r)[C(a , l (x-x i . 1 ),ri,^ a {xi-un) + 

g q( J »-*) 


(A. 75a) 


(A. 75b) 


+ e n S(<r*(x -Xi.i),r},fi)(f> g {xi,-T])^ + 

+ Mo e "“a ^“ =1 A ’C (<7‘(x - Xi-i), & p). 

To obtain the proper relations for an isotropic source, define the boundary 
conditions at the slab boundaries as: 


Leftmost slab - Isotropic source on left face 


F£ l (*0 = S 9 , 


(A. 76a) 


F 9 R{p) = 4>g{x l,-p). 


(A. 76b) 


Rightmost slab - Vacuum on right face 


Fl NS {p) = <t>g{x NS-li/0 + S%e i^*=i A », 


(A. 77a) 


f 3 r ’ ns (m) = 0. 


(A. 77b) 
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Middle slabs - Collided plus uncollided on both faces 


FlV) = <f> g (xi-i,n) + S$e 


F r'(m) = 


(A. 78a) 
(A. 78b) 


A pattern is set up when the boundary values are substituted into Rl^(/i,x) 
and R2j(/i,x) which is 

= J Q dl l 7 [C(<T*(Xj - x),Tl,fi)<j> g (Xi,-T)) + 


+ e ” S(<t*(x,- - x),r},n)<f> g {xi- U Tf) 


+ 


(A. 79a) 


+ Sq J^drj rje A 5 +<r i (x ri-l) )s((Tg(x t - - x),t),h), 

R- 2 ^’ 1 ) = Jo dr > * [ C K( X ~ + 

+ e " S(cr*_ 9 (x - x,_x), T), fi)<f> g {xi, -t}) 

+ So J Q dT] rje~^ £-*=i A ®C(cr‘(x - fi). 


(A. 79b) 


The Fn approximation could be used to put every flux in terms of F^ l (pt) 
and F^ NS (/i), but the result is a multiplicative recursive relationship or matrix 
over all the slabs. For any multiple slab analysis, it will be computationally faster 
to use the post processor to calculate the fluxes at the integration points used for 
RV g {n,x) and R2 ' g {n,x). 


A. 7 Semi-infinite Media 

This section is not a detailed derivation of the semi-infinite media Fn 
method, but the resultant equations are shown for completeness. 
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A. 7.1 Matrix Equations 


The boundary matrix equation is 
N - 1 1 

a a B ®(^) = Rl,(t',X 0 ) + Tl ff (l/,X 0 ,00) + 


Of=0 


<r a ~ 


fl — 9 


— —Slg(l/, x 0 ,oo). 

a 9-+g 

The interior matrix equations are 
N-l 


£ = Rl g (v,Xj) + — 1 T1 9 (i^, X j,oo)+ 

cr=0 a 9—9 


£ [d^’ J B^(i/) - c™A* (i,)] = R2 sl (t', Xj) H T2 3 (i/, x 0 ,Xj)+ 

a=0 


j “ — Sl ff (l/, Xj,00), 


t — — S2 a (i/, x 0 ,x ; ), 


9~*9 


where, 

Rl 3 (i/, Xj) = ^ 


Sq — — — e~ Za ^ x Xo \ BeamSource 

Ho + v 

Sq f dr/ — - — e V" (*-*o)^ IsotropicSource, 
l ./o T) + 1/ 


R2 s (j/,x i ) = < 


So/ioC(<r 3 (x — x 0 ), /io, ^), BeamSource 


| Sg J di] T)C(<7 g (x — Xo),Tj,i;), IsotropicSource, 


g u g 1/ 


C -i 


V ^ v 


-e * r) = v, 


9-1 


Tl 3 (i/, x j , oo ) = (T g 0 y_ s I* a (i/, x_j , oo), 


3 ' = 1 


5-1 


3\j) ^ ] ^3'— ’5^33 (^> 3Jj)i 

3 ' = 1 


(A. 80a) 

(A. 80b) 
(A. 80c) 

(A.80d) 

(A.80e) 

(A.80f) 

(A.80g) 


(A.80h) 
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t/ = 0 


Sl s (t/,Xj,oo) = 


Sj(Xj) 

^ /°°dz a (z) 1/^0, 


and 


S2 fl (z/,x 0 ,oo) = < 


A. 7. 2 Post Processor 


i/ = 0 


Fdz c -?(-i-)S f (z) *#0. 

J/ Jro 


The boundary post processor is 
1 


4>g (XO) /*) — 


2 *A(aO L 


I] a£BaV)+ 

or=0 


"h ^ Rig (jCx^xq) 4" T l^f (/^? xo> l '~ ) o ) 4“ S ^ Xoi )] * 
The interior post processor is 
1 




jV-1 


E {^BM +d«A>(M)} + 


a=0 


2^5 ^$0 (^) 

+cr^Rl^(/i, Xj) + Tlj(jz, Xj, oo) + SI g (fi, Xj, oo)] , 


<^ 5 (xj, /i) — 


iV-1 


E {dS'BJM + (cj-i- 


2*»*»00 [ „ 0 
_ a ».. e -?(*,-»>) A»(m)} + <r,-,R2,0i,ij)+ 

+T2 ff (/z, x 0 , Xj) + S2g(^,x 0 ,x J )] . 

A. 7.3 Scattering Terms 

The energy group particle transfer terms for £ 6 [0,s ffff /] are 

AT-l i 

X i 00 ) = Sl 5 '(s^/^, X, 00) + 


9‘~ 1 


a=0 


<y’g'^g ( 


b — 1] <V'- s ' I ss"(f> ;r > 00 )> 


(T g'~'9' g" = \ 


(A.80i) 


(A.80j) 


(A. 81a) 


(A. 81b) 


(A. 81c) 


(A. 82a) 
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n-i / 1 

°°) = < £ J ^«( s w , 0 - S2 A s gA> x > °°)+ 

a=0 " ' 

g'-i 


Gg'->9* 


)•— ®- E "V-WJJ-tt. *.<»)> 


a a'-*9' g" = l 


where, 


K = < 


a" x = x 0 
c£’ J x ^ x 0 . 


The energy group particle transfer terms for £ ^ [0, s 33 <] are 


N - 1 




a=0 


, ni / c ~\ ( 3 99 t €)I J { 3 93 , 0 Q^ / £ _ \ 




(T gt 


a g X^ , £, oo), 


s"=l 


iV-1 


ff.A,v(»«'f)JS'«.*o.*) = E [d" A* (-*„.{) + (c“- 


a=0 


-a^e s » ,<(r Xo) ) A 9 a {s ggl {) 


+ 


i DO /« £ _\ ( 6 P5 / 0L(-5 55 ^) Co ( ^ 

+ , x) ? ^o? £ j 

CTg> 

g f - 1 

( 5 0S'O ^sV'^ 7 ^Oi ^ )j 

<7" = 1 


where, 

ACx = { 


3- a X — Xq 

cl’> X £ x 0 , 


/Co = \ 


0 x = xq 


(A. 82b) 


(A. 82c) 


(A. 83a) 


(A. 83b) 


(A. 83c) 


d£ J x ^ Xq. 


(A. 83d) 
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A. 7. 4 Scalar Flux 


The scalar flux is 
ff-1 


M x <i) = 

*,(*) = ^ E (<£* + <•’) Mi) + 1, 


where, 


I = l 


-ZL(x- 


Soe •‘o 

l *£« 


(x-x 0 ) 


BeamSource 


e x ° * IsotropicSource. 


A. 7. 5 GlT Neutron Source 

The neutron source from the galactic cosmic ray cascade is 
„ J-2 Nt g - 


S gi x ) = Y, vfj 


^O/iGIT — x “i l, ‘ g^x S />( E <?') 

where, 

T / \ f°° A — (r— u/) — (<7 j_, — y_ t1 )u/— t^-s/ 

/(x) = — awe ^git v ; v J ’i b r 7 11 ; ^ } ; 

J Xn 


I(x), 


' Xq 

S a i 


\ VJ-ii Vj-i 2 / 


SI 


The galactic cosmic ray cascade neutron source in the Fn context i 

J — 2 N t 

MqJ-iCTj-i 

^CToMGIT /=1 
J—2 




r=l 

N, 


e,s P (E,.) 


S2 a (p,x„,x) = 

^O/iGIT ^ Sp(E 3 .) 


/2(x), 


where, 


71(x) = j™&ze- a -t {z ~ x) I(z), 


(A. 84a) 
(A. 84b) 

(A. 84c) 

(A. 85a) 

(A. 85b) 

(A. 86a) 
(A. 86b) 


(A. 86c) 
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I2(x) = f X dz e~^ {x ~ z) I(z). (A.86d) 

Jx o 

The integrals are evaluated analytically using MathCAD. 

A. 8 Critical Slab Width and Associated Flux 
From the Fn formulation, a critical slab width and flux can be determined 
for the one slab, one energy group, source free problem. The critical problem is 
defined by the infinite medium dispersion relation with an imaginary root and the 
specification of the resultant matrix from collocation being equal to zero. 
Normalization of the angular flux is an input parameter, and the flux is considered 
to be symmetric about the center of the slab. 

A. 8.1 Critical Slab Width 

The boundary matrix equations with one group, one slab, no source, and 
flux symmetry ( a a = 6 a ) reduce to 

E + A a (i/)e~^] = 0. (A. 87) 

a=0 

To determine j/o, use the imaginary root of the infinite medium dispersion relation 

A(i'o) = 0 = l — i/o — tan -1 ( — ^ . (A. 88) 

a \i/oJ 

Since i/q is purely imaginary, the first row of the expansion coefficient matrix is 
complex. To find the critical width, the real and imaginary parts of the 
determinant of the matrix must be zero at the same value of A c (the critical 
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width). With /3 = 1, . . . , N — 1 and a — 0, . . . , N — 1, the determinant equations 
are 


3ftct0 

B a {vp) + A a (up)e u » 


= 0, 


(A. 89a) 


B a {vp) + A a (up)e u e 


= 0 , 


(A. 89b) 


where, 


3£ Q o = SA a (cos — — l) + z 0 SB o sin — -1 , 

la L V z 0 / zq J 

SA q sin — + z 0 SB a ( cos — + 1^1 , 

L ZQ \ Zq J J 


o 

= 2a 


(A. 89c) 
(A.89d) 


and, 


i/ 0 = iz 0 , (A.89e) 

A c = a (x c — 0) , (A.89f) 

SA a = / drj 2 ^ a (??), (A.89g) 

SB a = f dr/ - Y 2 0 a (7?). (A.S9h) 

./o TJ + ^0 


The determinant is found from the LU decomposition routine in Reference [24]. A 


bisection method is used to find A c to satisfy the above relations. 
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A. 8. 2 Boundary and Interior Critical Flux Values 

To determine the boundary flux, rewrite the matrix equation as 

+ A a (u 0 )e~^} = -a o [Bo(iy 0 ) + A o (u p )e- A ^}, (A.90) 

0 = 1 

where, ao is a normalization parameter. To calculate the additional a a values, 
invert the matrix and use the normal post processor to find the fluxes. 

To find the interior flux, the matrix equations for (3 = 0, 1, . . . , /V — 1 are 


N-l 

T [c a B a M-d a A a (v p )} = 

N-l 

_ e -<r{x c -x)/v 0 a a A a (v 0 ), 

(A. 91a) 

ot=0 

a=0 



N - 1 


[daB a (vp) - C a A a {ug)] = 

Of = 0 

_ e -*(*-0 )/u p £ a a A a (i/p). 

a=0 

(A. 91b) 


However, for /3 = 0, the A a (i/pY s and B a (vpYs are complex 


A a (i/p) = — 1 [SA* - iz 0 SB Q ] , 

Z(T 


B a {vp) = ^ [SA a + i^oSBa] 

Id 


(A. 92a) 
(A. 92b) 


The expansion coefficients are real, but the matrix terms and the matrix 
inversion routine need to be complex in order to calculate the expansion 
coefficients properly. Once the c a and d a values are found, the normal post 
processor can be used to determine the interior flux distribution. 


A. 8. 3 Fn to Anisn/Pc Normalization 

To calculate the normalization factor between the Fn method and 
ANISN/Pc, set the scalar flux found at x = 0 from both methods equal. The 
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scalar flux at the boundary for the Fn method is 

<M0) = p-°o- (A. 93) 

Zcr 

For a certain set of cross sections, vcrj + a a = cct, set the normalization factor used 
in the Fn algorithm as 

ScaIe = 2< ftanisn(0) (A. 94) 

c 
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Appendix B 

Program User’s Manual 


This chapter provides a user’s manual for the operation of MGSLAB, 
MGSEMI, and FNCRIT. The input to and the output from the programs will be 
discussed. An example is included for clarity. The programs can be acquired 
through the University of Arizona’s Nuclear Engineering Department. 

Four versions of the program were created in the course of the research on 
the Fn method. The first version analyzed finite slab boundaries only, but 
included a degenerate eigenvalue algorithm. This version exists only to verify the 
algorithm and final program against existing scientific literature (References [11] 
and [12]). A second version of the program, referred to as MGSLAB, analyzes 
boundaries and interior slab positions without the degenerate eigenvalue algorithm 
which usually does not appear when normal material cross sections are used. This 
is the workhorse of the four programs. An algorithm that analyzes one group, 
homogeneous critical slabs is the third program created. It is referred to as 
FNCRIT. The last version is a homogeneous semi-infinite variant called MGSEMI. 
These programs use common subprograms shared in an object library. 
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B.l Program Input 

The input to each program is the same, whether or not the values are used 
by the particular program version being executed. The only exception is the 
scaling value used by FNCRIT. The pseudo-code below outlines the input flow. 

The programs read this file from FORTRAN unit 5. Under most operating systems, 
this unit can be defined as a file name outside the program. 

Line 1: 

Nslabs - Number of slabs [negative for diagnostic] 
(<MAXSLAB). 

Ngroup - Number of energy groups (<MAXGROUP). 

Inner - Number of boundary inner iterations for each slab at 
each FN iterate. 

GLBoun - Inhomogeneous term integration quadrature 
(<MAXQUAD) 

Line 2: 

Nstart - Initial N for FN iteration: 

• Negative value for collocation points of Chebyshev 
roots. 

• Positive value for collocation points of Legendre 
roots. 

Nend - Final N for FN iteration (<MAXN) 

Nstep - Increment between approximations 

Line 3: 

EditMu - Number of positive direction edit points 
(<MAXEDIT) 

EMType - Type of direction edit points to use: 

1. Gauss-Legendre points of order EditMu 

2. Evenly spaced points 
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3. User supplied point values 

If EMType .eq. 1 then Line 4: 

EMstrt - Starting edit direction value (0.0) 

EMend - Final edit direction value (1.0) 

If EMType .eq. 2 then Line 4: 

EMstrt - Starting edit direction value (0.0) 

EMend - Final edit direction value (1.0) 

EMZero - The value to use for a Mu of zero (1.0d-10) 

If EMType .eq. 3 then Line 4: 

EMuser - User input positive edit directions (L.EditMu) 

Line 5: 

Tol - Flux convergence tolerance (1.0d-5) 

Line 6: 

Stype - The type of external source 

0. No source 

1. Beam source at left face of slab 

2. Isotropic source at left face of slab 

If Stype .eq. 1 then Line 7: 

MuO,SO - Direction and Intensity of external source for each 
group (L.Ngroup) 

If Stype .eq. 2 then Line 7: 

SO - Intensity of external source for each group (L.Ngroup) 

Line 8: 


Gsour - Number of energy groups that have distributed 
sources (if a -1, use the analytical GIT source) 


If Gsour .It. 0 then Line 9: 
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MuG - Direction of incoming GIT Beam 
dens - Density of the target material 


Line 10: 


Width - Slab width in centimeter (ignored for the FNCRIT 
program). 

EditX - Number of spatial edit points in this slab 
(<MAXEDIT) 

EXType - Type of spatial edit points used 

1. Evenly spaced points within slab (<MAXEDIT) 

2. User supplied points within slab (< MAXED IT) 

3. Gauss-Legendre points of order EditX within slab 
(<MAXEDIT) 

GLMatx - Gauss-Legendre integration order for evaluation 
of integrals in the matrix terms for this slab 
(<MAXQUAD) 

BType - Type of basis functions to use in this slab 

1. Shifted Legendre: P/( 2x — 1) 

2. Modified Shifted Legendre: P/(2(x md — 1) 

3. Monomials: x* 

4. Modified Shifted Monomials: (2x md — 1)* 

md - Parameter used for BType 2 and 4 (0.75) 

Line 11: 

Sg,Sgpg - Total and down scatter cross sections. Array 

dimensions for down scatter cross sections: (slab, from 
group, to group). 

st ( 1 ) sl(l,l) 0 0 ... 

st (2) s2 ( 1 , 2) s2(2,2) 0 ... 

st (3) s3( 1 , 3) s3(2 ,3) s3(3,3) . . . 


If EXType .eq. 2 then Line 12: 


EXuser - User input spatial edit points (1.. EditX) 
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Line 13 (used by FNCRIT only): 

Scale - The flux scaling value 
Last Line optional: 

PNodes - Optional nodes to include in the boundary flux 
plot file (<100) 

If the GlT source is used, two more file must be read. FORTRAN unit 40 
contains the Yf“ r ( values generated by the GlT program. These are in binary 
format. FORTRAN unit 41 contains the input deck used by the GlT program. 

Line 1: 

JJ - Maximum species (charge) number (<MAXJ) 

LJJ - Charge number of last species 

Line 2: 

IRM - 

-1. Use wilson stopping power 
1. Use simplified stopping power 
TOL - Relative error for romberg numerical integration 

Line 3: 

IDFL - 1 For flux profile for incident beam of species JJ 

(Only option used) 

ISX - 

-1. Use edit flux versus x 
1. Use edit flux versus path length 

Line 4: 

IB IN - 

-1. Read binary tape 40 for y coefficients 
0. No action 
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1. Write binary tape 40 for y coefficients 

IPR - 

-1. All output on one file (unit 21) 

1. Output on multiple files (units 21 through 39) 

Line 5: 

EO - Input beam initial energy 

LE - Number of energy edit points (<20) 

LX - Number of spatial edit intervals (< 100) 

Line 6: 

SO - Model cross section normalization 
FN - Fraction of neutrons produced in an interaction, /„ 
Ell - First energy for simplified stopping power 
E22 - Second energy for simplified stopping power 

If ISX .It. 0 then Line 7: Energy edit grid (EN(I),I=1,LE) 

If ISX .gt. 0 then Line 7: Spatial grid (X(I),I=1,LX) 

Section B.3 shows an example of an input deck. 


B.2 Program Output 

The output routines generate various types of output data. The list below 
shows FORTRAN unit number written and the type of data written to it. Under 


most operating systems, these units can be defined as file names outside the 


program. 

• Unit 6: the Screen 

• Unit 7: diagnostics file 

• Unit 21: node and group angular fluxes in text format 
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• Unit 22: slab and group angular fluxes in text format 

• Unit 23: group scalar fluxes in text format 

• Unit 24: total scalar flux in text format 

• Unit 25: total scalar flux in plot format 

• Unit 26: screen echo to a file 

• Unit 27: node and group angular fluxes in plot format 

• Unit 28: slab and group angular fluxes in plot format 

• Unit 29: group scalar flux in plot format 

• Unit 30: the GlT ion flux in plot format 

• Unit 31: the GlT neutron source in plot format 

• Unit 49: energy group, slab position, and group scalar flux in three 
dimensional plot format 

• Unit 50 to 69: energy group, direction, and group angular flux in three 
dimensional plot format 

Section B.3 shows the text output from the example input deck. 

B.3 Example 

This section shows a typical example input deck and the associated output 
text files. The physical situation being modeled is a three energy group, three slab 
problem. A beam source of unit strength for group one only is incident on the left 


face of the media. 



Input Deck: 


3 3 3 35 
11 49 2 
16 1 
0.0 1.0 
1.0d-04 
1 

1.0 1.0 1.0 0.0 1.0 0.0 
0 

0.2d0 9 1 35 2 0.75d0 
2 . OdO l.OdO O.OdO O.OdO 
5 . OdO 1 .OdO 3. OdO O.OdO 
3. OdO O.ldO 0 . 2d0 1.5d0 
0.4d0 19 1 35 2 0.75d0 
5. OdO 3. OdO O.OdO O.OdO 
2. OdO 0 . 5d0 l.OdO O.OdO 
3. OdO 0.4d0 O.ldO l.OdO 
0 . 4dO 19 1 35 2 0.75d0 
3. OdO 2 . 8d0 O.OdO O.OdO 
2. OdO O.ldO 0 . 9d0 O.OdO 
5. OdO 0 . 9d0 0 . 9d0 l.OdO 
5 20 40 


\ Nslabs Ngroup Inner GLBoun 
\ Nstart Nend Nstep 
\ EditMu EMType 
\ EMstart EMend 
\ Tol 

\ Beam Source 

\ MuO SO for groups 1 , 2 , and 3 
\ No GIT Source 

\ Width EditX EXType GLMatx BType md 
\ Cross sections 

\ Width EditX EXType GLMatx BType md 
\ Cross sections 

\ Width EditX EXType GLMatx BType md 
\ Cross sections 

\ PNodes 
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Unit 26 - screen echo to a file: 

3D Scalar Plot Data File Generated 
Energy Group * 0001 

N Conv ToCon Pos max err Neg max err Pos ave err Neg ave err 

0011 0000 0096 1 . 000000E+00 1.000000E+00 1.000000E+00 1.000000E+00 

0011 0000 0096 1 . OOOOOOE+OO 1.000000E+00 1.000000E+00 1.000000E+00 

0011 0000 0096 1. 000000E+00 1. 000000E+00 1. 000000E+00 1. 000000E+00 

0013 0000 0096 3.553633E-03 5.712922E-03 1.167111E-03 1.301876E-03 

0013 0000 0096 3.872581E-03 6.227630E-03 1.272876E-03 1.419210E-03 

0013 0000 0096 3.901608E-03 6.274470E-03 1.282503E-03 1.429885E-03 

0015 0093 0096 1.489162E-04 1.064508E-04 1.206050E-05 6.633653E-06 

0015 0093 0096 1.488493E-04 1.064160E-04 1.209246E-05 6.696274E-06 

0015 0093 0096 1.488428E-04 1.064128E-04 1.209542E-05 6.701626E-06 

0017 0096 0096 3.640764E-05 2.877266E-05 4.588094E-06 2.083833E-06 

Boundary Flux Converged 

0017 0096 1600 1. 000000E+00 1. OOOOOOE+OO 9.400003E-01 9.400001E-01 

0019 1584 1600 1.766174E-04 1.468803E-04 6.441155E-06 4.480234E-06 

0021 1596 1600 1 . 158229E-04 1.116657E-04 1.090374E-06 9.166960E-07 

0023 1600 1600 2.720663E-05 3.480179E-05 1.227406E-07 1.803098E-07 

All Fluxes Converged 

Including 0.100000000000000 in the boundary flux output 

Including 0.400000000000000 in the boundary flux output 

Including 0.800000000000000 in the boundary flux output 

Energy Group = 0002 

N Conv ToCon Pos max err Neg max err Pos ave err Neg ave err 

0011 0000 0096 1. 000000E+00 1. OOOOOOE+OO 1. 000000E+00 1. 000000E+00 

0011 0000 0096 1. 000000E+00 1. OOOOOOE+OO 1. OOOOOOE+OO 1. 000000E+00 

0011 0000 0096 1. 000000E+00 1. 000000E+00 1. 000000E+00 1. 000000E+00 

0013 0000 0096 9.513772E-04 1.639784E-03 4.420175E-04 4.710332E-04 

0013 0000 0096 1.016245E-03 1.755531E-03 4.723103E-04 5.034820E-04 

0013 0000 0096 1.020552E-03 1.763143E-03 4.743243E-04 5.056349E-04 

0015 0094 0096 1.350178E-04 2.991238E-05 1.619709E-05 3.357026E-06 

0015 0094 0096 1.350328E-04 2.992016E-05 1.619631E-05 3.368438E-06 

0015 0094 0096 1.350335E-04 2.992069E-05 1.619635E-05 3.369010E-06 

0017 0096 0096 1.252745E-05 3.572194E-06 1.847479E-06 5.871654E-07 

Boundary Flux Converged 

0017 0096 1600 1. 000000E+00 1. 000000E+00 9.400001E-01 9.400000E-01 

0019 1590 1600 1.325581E-04 1.302121E-04 5.260669E-06 4.477305E-06 

0021 1600 1600 5 . 069402E-05 4.537826E-05 3.722857E-07 5.015554E-07 

All Fluxes Converged 

Including 0 . 100000000000000 in the boundary flux output 
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Including 0.400000000000000 in the 
Including 0.800000000000000 in the 
Energy Group * 0003 


Conv ToCon Pos max err 


N 

0011 0000 0096 
0011 0000 0096 
0011 0000 0096 
0013 0063 0096 
0013 0063 0096 
0013 0063 0096 
0015 0096 0096 


1 . 000000E+00 
1 . OOOOOOE+OO 
1 . 000000E+00 
7 . 529517E-04 
7 . 670185E-04 
7.673035E-04 
1 . 967 129E-05 


Neg max err 
1. 000000E+00 
1. OOOOOOE+OO 
1. 000000E+00 
1 . 363943E-03 
1.400001E-03 
1.400694E-03 
2 . 913212E-05 


Boundary Flux Converged 


0015 0096 1600 1. 000000E+00 1. 000000E+00 
0017 1590 1600 1 . 684172E-04 1.991901E-04 
0019 1600 1600 5.225291E-05 6.332041E-05 


All Fluxes Converged 

Including 0.100000000000000 in the 
Including 0.400000000000000 in the 
Including 0.800000000000000 in the 


boundary flux output 
boundary flux output 


Pos ave err 
1. 000000E+00 
1. OOOOOOE+OO 
1. 000000E+00 
1 . 882021E-04 
1.922363E-04 
1 . 923174E-04 
2.491874E-06 

9.400001E-01 
4.851583E-06 
6 . 616150E-07 


Neg ave err 
1. 000000E+00 
1. 000000E+00 
1. 000000E+00 
2 . 965920E-04 
3 . 036932E-04 
3 . 038329E-04 
1 . 872452E-06 

9.400001E-01 
3.744157E-06 
5 . 106637E-07 


boundary flux output 
boundary flux output 
boundary flux output 
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Unit 21 - node and group angular fluxes in 
some of the spatial positions) 

text format: (only including 

Angular flux at the Direction 

Edit Points 

for group 1 

Spatial point 0.00000 

Mu 

Mu < 0 

Mu > 0 

5 . 29953E-03 

3 . 219E-01 

0 . 000E+00 

2.77125E-02 

3. 218E-01 

0 . 000E+00 

6.71844E-02 

3 . 188E-01 

0 . 000E+00 

1 . 22298E-01 

3 . 146E-01 

0 . 000E+00 

1 . 91062E-01 

3 . 103E-01 

0 . 000E+00 

2 . 70992E-01 

3.048E-01 

0 . 000E+00 

3.59198E-01 

2 . 977E-01 

0 . 000E+00 

4 . 52494E-01 

2.895E-01 

0 . 000E+00 

5 . 47506E-01 

2 . 808E-01 

0 . 000E+00 

6.40802E-01 

2.724E-01 

0 . OOOE+OO 

7 . 29008E-01 

2.647E-01 

0 . 000E+00 

8 . 08938E-01 

2 . 579E-01 

0. 000E+00 

8 . 77702E-01 

2 . 522E-01 

0. 000E+00 

9 . 32816E-01 

2.478E-01 

0. 000E+00 

9 . 72288E-01 

2.447E-01 

0. 000E+00 

9 . 94700E-01 

2.429E-01 

0. 000E+00 

Spatial point 0.02000 

Mu 

Mu < 0 

Mu > 0 

5 . 29953E-03 

3 . 219E-01 

3 . 222E-01 

2 . 77125E-02 

3 . 201E-01 

2.463E-01 

6.71844E-02 

3 . 159E-01 

1.447E-01 

1 . 22298E-01 

3. 115E-01 

8 . 995E-02 

1 . 91062E-01 

3 . 074E-01 

6 . 091E-02 

2 . 70992E-01 

3.020E-01 

4.425E-02 

3 . 59198E-01 

2 . 948E-01 

3 . 398E-02 

4 . 52494E-01 

2 . 864E-01 

2 . 728E-02 

5 . 47506E-01 

2.777E-01 

2 . 272E-02 

6.40802E-01 

2 . 692E-01 

1.951E-02 

7 . 29008E-01 

2 . 614E-01 

1.721E-02 

8 . 08938E-01 

2.546E-01 

1 . 556E-02 

8 . 77702E-01 

2.489E-01 

1.436E-02 

9 . 32816E-01 

2.445E-01 

1.353E-02 

9.72288E-01 

2.414E-01 

1.300E-02 

9 . 94700E-01 

2 . 397E-01 

1.271E-02 


(edited file here) 
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Spatial point 0.98000 
Mu 

5 . 29953E-03 
2.77125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2.70992E-01 
3 . 59198E-01 
4.52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 
Spatial point 1.00000 
Mu 

5 . 29953E-03 
2 . 77125E-02 
6 . 7 1844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5 . 47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9.32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Angular flux at the Direction 
Spatial point 0.00000 
Mu 

5 . 29953E-03 
2 . 77 125E-02 


Mu < 0 

Mu > 0 

6 . 656E-02 

6.792E-02 

5.714E-02 

7.062E-02 

3.763E-02 

7.505E-02 

2.460E-02 

8 . 084E-02 

1 . 706E-02 

8 . 744E-02 

1 . 256E-02 

9 . 396E-02 

9 . 723E-03 

9 . 978E-02 

7.846E-03 

1 . 049E-01 

6 . 556E-03 

1 . 095E-01 

5 . 644E-03 

1 . 134E-01 

4 . 989E-03 

1.167E-01 

4 . 513E-03 

1.193E-01 

4.171E-03 

1 . 212E-01 

3 . 933E-03 

1 . 226E-01 

3.778E-03 

1 . 234E-01 

3 . 695E-03 

1 . 239E-01 

Mu < 0 

Mu > 0 

0 . 000E+00 

5 . 925E-02 

0 . 000E+00 

6 . 262E-02 

0 . 000E+00 

6 . 760E-02 

0 . 000E+00 

7 . 381E-02 

0. 000E+00 

8 . 081E-02 

0 . 000E+00 

8 . 779E-02 

0 . OOOE+OO 

9 . 412E-02 

0 . OOOE+OO 

9 . 970E-02 

0. 000E+00 

1 .046E-01 

0. 000E+00 

1 . 089E-01 

0 . 000E+00 

1 . 125E-01 

0. 000E+00 

1 . 153E-01 

0. OOOE+OO 

1 . 174E-01 

0 . 000E+00 

1 . 189E-01 

0 . 000E+00 

1 . 198E-01 

0 . OOOE+OO 

1 . 203E-01 

Edit Points for 

group 2 

Mu < 0 

Mu > 0 

1 . 886E-01 

0. 000E+00 

1 . 932E-01 

0. 000E+00 
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6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2.70992E-01 
3 . 59198E-01 
4.52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Spatial point 0.02000 
Mu 

5 . 29953E-03 
2.77125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5.47506E-01 
6 . 40802E-01 
7 . 29008E-01 
8.08938E-01 
8.77702E-01 
9 . 32816E-01 
9.72288E-01 
9 . 94700E-01 

(edited file here) 

Spatial point 0.98000 
Mu 

5 . 29953E-03 
2.77125E-02 
6 . 7 1844E-02 
1 . 22298E-01 
1 . 91062E-01 


1 . 985E-01 

0 . 000E+00 

2 . 032E-01 

0 . 000E+00 

2 . 067E-01 

0.000E+00 

2 . 081E-01 

0 . 000E+00 

2 . 070E-01 

0 . OOOE+OO 

2 . 035E-01 

0 . 000E+00 

1 . 983E-01 

0. 000E+00 

1 . 924E-01 

0. 000E+00 

1 . 864E-01 

0. 000E+00 

1 . 808E-01 

0. 000E+00 

1 . 760E-01 

0. 000E+00 

1 . 722E-01 

0. 000E+00 

1 . 696E-01 

0. 000E+00 

1 . 681E-01 

0. 000E+00 

Mu < 0 

Mu > 0 

2 . 050E-01 

2 . 000E-01 

2 . 07 IE-01 

1 . 969E-01 

2 . 098E-01 

1 . 541E-01 

2 . 123E-01 

1 . 108E-01 

2 . 137E-01 

8 . 066E-02 

2 . 132E-01 

6 . 102E-02 

2 . 102E-01 

4 . 802E-02 

2 . 051E-01 

3 . 917E-02 

1.986E-01 

3 . 297E-02 

1 . 916E-01 

2.853E-02 

1.848E-01 

2 . 531E-02 

1.786E-01 

2 . 296E-02 

1.735E-01 

2 . 126E-02 

1 . 694E-01 

2 . 007E-02 

1 . 666E-01 

1 . 929E-02 

1 . 650E-01 

1 . 888E-02 


Mu < 0 Mu > 0 

1.425E-02 1.458E-02 

1 . 062E-02 1 . 526E-02 

6 . 194E-03 1 . 644E-02 

3 . 843E-03 1 . 822E-02 

2 . 599E-03 2 . 141E-02 



217 


2.70992E-01 
3 . 59198E-01 
4. 52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 
9.32816E-01 
9.72288E-01 
9 . 94700E-01 

Spatial point 1.00000 
Mu 

5 . 29953E-03 
2 . 77 125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9.32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Angular flux at the Direction 

Spatial point 0.00000 
Mu 

5 . 29953E-03 
2 . 77 125E-02 
6 . 7 1844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5.47506E-01 
6.40802E-01 


1 . 887E-03 

2 . 700E-02 

1.449E-03 

3.481E-02 

1.163E-03 

4.358E-02 

9 . 684E-04 

5. 194E-02 

8.317E-04 

5.905E-02 

7 . 337E-04 

6 . 463E-02 

6 . 630E-04 

6.876E-02 

6 . 122E-04 

7 . 166E-02 

5 . 768E-04 

7 . 360E-02 

5 . 538E-04 

7.479E-02 

5.416E-04 

7 . 540E-02 

Mu < 0 

Mu > 0 

0 . 000E+00 

1.320E-02 

0 . 000E+00 

1 . 397E-02 

0 . 000E+00 

1 . 519E-02 

0 . 000E+00 

1 . 696E-02 

0 . 000E+00 

1.995E-02 

0 . 000E+00 

2 . 518E-02 

0 . 000E+00 

3 . 259E-02 

0 . 000E+00 

4. 105E-02 

0 . OOOE+OO 

4 . 925E-02 

0 . 000E+00 

5 . 631E-02 

0. OOOE+OO 

6 . 192E-02 

0. OOOE+OO 

6.610E-02 

0. OOOE+OO 

6 . 908E-02 

0. 000E+00 

7 . 108E-02 

0. 000E+00 

7 . 233E-02 

0. OOOE+OO 

7 . 297E-02 

Edit Points 

for group 3 

Mu < 0 

Mu >0 

4 . 235E-02 

0. 000E+00 

4.407E-02 

0. OOOE+OO 

4 . 606E-02 

0. OOOE+OO 

4 . 823E-02 

0. OOOE+OO 

5 . 096E-02 

0. OOOE+OO 

5.385E-02 

0. OOOE+OO 

5 . 619E-02 

0. OOOE+OO 

5 . 773E-02 

0. OOOE+OO 

5 . 853E-02 

0. OOOE+OO 

5.878E-02 

0. OOOE+OO 



7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 
9.32816E-01 
9.72288E-01 
9 . 94700E-01 

Spatial point 0.02000 
Mu 

5 . 29953E-03 
2.77125E-02 
6 . 71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4.52494E-01 
5 . 47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 


5.866E-02 

0 . 000E+00 

5 . 834E-02 

0 . 000E+00 

5.794E-02 

0 . 000E+00 

5.756E-02 

0 . 000E+00 

5.725E-02 

0 . 000E+00 

5.707E-02 

0 . 000E+00 

Mu < 0 

Mu > 0 

4.661E-02 

4. 608E-02 

4.761E-02 

3 . 990E-02 

4.898E-02 

2 . 639E-02 

5 . 078E-02 

1.727E-02 

5 . 343E-02 

1 . 199E-02 

5 . 622E-02 

8 . 828E-03 

5.835E-02 

6 . 834E-03 

5 . 963E-02 

5.516E-03 

6 . 017E-02 

4 . 609E-03 

6 . 020E-02 

3 . 969E-03 

5 . 989E-02 

3 . 508E-03 

5 . 942E-02 

3.174E-03 

5 . 891E-02 

2.933E-03 

5 . 844E-02 

2.765E-03 

5 . 808E-02 

2 . 656E-03 

5 . 786E-02 

2 . 598E-03 


(edited file here) 


Spatial point 0.98000 
Mu 

5.29953E-03 
2.77125E-02 
6 . 7 1844E-02 
1 . 22298E-01 
1 . 91062E-01 
2.70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 


Mu < 0 

Mu > 0 

2 . 054E-02 

2 . 093E-02 

1 . 966E-02 

2 . 147E-02 

1 . 535E-02 

2 . 233E-02 

1 . 099E-02 

2.348E-02 

7 . 993E-03 

2 . 485E-02 

6 . 042E-03 

2.639E-02 

4 . 753E-03 

2 . 800E-02 

3 . 875E-03 

2 . 960E-02 

3 . 261E-03 

3. 111E-02 

2 . 822E-03 

3 . 244E-02 

2 . 503E-03 

3 . 354E-02 

2 . 270E-03 

3.440E-02 

2 . 102E-03 

3 . 502E-02 



9 . 32816E-01 
9.72288E-01 
9 . 94700E-01 

Spatial point 1.00000 
Mu 

5 . 29953E-03 
2.77125E-02 
6 . 7 1844E-02 
1 . 22298E-01 
1 . 91062E-01 
2.70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 


1 . 984E-03 

3 . 546E-02 

1 . 908E-03 

3 . 573E-02 

1 . 867E-03 

3 . 587E-02 

Mu < 0 

Mu > 0 

0 . 000E+00 

1 . 815E-02 

0 . 000E+00 

1.885E-02 

0 . 000E+00 

1 . 987E-02 

0 . 000E+00 

2 . 115E-02 

0 . 000E+00 

2.262E-02 

0 . 000E+00 

2 . 423E-02 

0 . OOOE+OO 

2 . 592E-02 

0 . OOOE+OO 

2.759E-02 

0. OOOE+OO 

2 . 916E-02 

0. 000E+00 

3 . 057E-02 

0. OOOE+OO 

3.174E-02 

0. OOOE+OO 

3 . 266E-02 

0. OOOE+OO 

3.335E-02 

0. OOOE+OO 

3.383E-02 

0. OOOE+OO 

3 . 414E-02 

0. OOOE+OO 

3.430E-02 



Unit 22 - slab and group angular fluxes in text format: 


Angular flux at the Direction Edit Points for group 1 


Spatial point 0.00000 
Mu 

5 . 29953E-03 
2.77125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5 . 47506E-01 
6 . 40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 
9 . 32816E-01 
9.72288E-01 
9 . 94700E-01 

Spatial point 0.20000 
Mu 

5 . 29953E-03 
2 . 77 125E-02 
6.71844E-02 
1 . 22298E-01 
1.91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5.47506E-01 
6 . 40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9.32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Spatial point 0.60000 
Mu 

5 . 29953E-03 


Mu < 0 

Mu >0 

3 . 219E-01 

0 . 000E+00 

3 . 218E-01 

0 . 000E+00 

3 . 188E-01 

0 . 000E+00 

3 . 146E-01 

0 . 000E+00 

3 . 103E-01 

0 . 000E+00 

3 . 048E-01 

0 . 000E+00 

2 . 977E-01 

0 . 000E+00 

2 . 895E-01 

0 . 000E+00 

2 . 808E-01 

0 . 000E+00 

2.724E-01 

0 . 000E+00 

2 . 647E-01 

0 . 000E+00 

2 . 579E-01 

0 . 000E+00 

2 . 522E-01 

0 . 000E+00 

2 . 478E-01 

0 . 000E+00 

2.447E-01 

0 . 000E+00 

2.429E-01 

0 . 000E+00 

Mu < 0 

Mu > 0 

3 . 331E-01 

2 . 786E-01 

3 . 298E-01 

2 . 813E-01 

3.233E-01 

2 . 856E-01 

3.140E-01 

2 . 806E-01 

3 . 025E-01 

2 . 589E-01 

2 . 897E-01 

2 . 296E-01 

2 . 768E-01 

2 . 008E-01 

2 . 645E-01 

1 . 760E-01 

2 . 532E-01 

1 . 558E-01 

2.433E-01 

1 . 397E-01 

2 . 346E-01 

1 . 272E-01 

2 . 274E-01 

1.176E-01 

2 . 215E-01 

1 . 104E-01 

2. 170E-01 

1 . 052E-01 

2 . 139E-01 

1 . 018E-01 

2 . 121E-01 

9 . 995E-02 

Mu < 0 

Mu > 0 

1 . 607E-01 

1 . 030E-01 
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2.77125E-02 

1 . 606E-01 

1 . 031E-01 

6 . 71844E-02 

1 . 592E-01 

1 . 041E-01 

1.22298E-01 

1.562E-01 

1.065E-01 

1 . 91062E-01 

1 . 516E-01 

1 . 103E-01 

2.70992E-01 

1.452E-01 

1 . 160E-01 

3 . 59198E-01 

1 . 374E-01 

1 . 229E-01 

4.52494E-01 

1 . 290E-01 

1 . 301E-01 

5.47506E-01 

1 . 207E-01 

1 . 364E-01 

6.40802E-01 

1 . 132E-01 

1.412E-01 

7 . 29008E-01 

1 . 067E-01 

1 . 445E-01 

8 . 08938E-01 

1 . 012E-01 

1 . 465E-01 

8 . 77702E-01 

9 . 693E-02 

1.476E-01 

9 . 32816E-01 

9 . 369E-02 

1 . 481E-01 

9 . 72288E-01 

9 . 149E-02 

1.482E-01 

9 . 94700E-01 

9 . 027E-02 

1.482E-01 

Spatial point 1.00000 



Mu 

Mu < 0 

Mu > 0 

5 . 29953E-03 

0 . 000E+00 

5 . 925E-02 

2.77125E-02 

0 . 000E+00 

6 . 262E-02 

6.71844E-02 

0 . 000E+00 

6.760E-02 

1 . 22298E-01 

0 . OOOE+OO 

7 . 381E-02 

1 . 91062E-01 

0 . 000E+00 

8 . 081E-02 

2 . 70992E-01 

0. 000E+00 

8.779E-02 

3 . 59198E-01 

0. 000E+00 

9.412E-02 

4 . 52494E-01 

0. OOOE+OO 

9 . 970E-02 

5.47506E-01 

0. OOOE+OO 

1 . 046E-01 

6.40802E-01 

0. OOOE+OO 

1 . 089E-01 

7 . 29008E-01 

0. OOOE+OO 

1 . 125E-01 

8 . 08938E-01 

0. OOOE+OO 

1 . 153E-01 

8 . 77702E-01 

0. OOOE+OO 

1 . 174E-01 

9 . 32816E-01 

0. OOOE+OO 

1 . 189E-01 

9 . 72288E-01 

0. OOOE+OO 

1 . 198E-01 

9 . 94700E-01 

0. OOOE+OO 

1 . 203E-01 

.ar flux at the Direction 

Edit Points for 

group 2 

Spatial point 0.00000 



Mu 

Mu < 0 

Mu > 0 

5 . 29953E-03 

1 . 886E-01 

0. OOOE+OO 

2 . 77 125E-02 

1 . 932E-01 

0. OOOE+OO 

6 . 7 1844E-02 

1 . 985E-01 

0. OOOE+OO 

1.22298E-01 

2 . 032E-01 

0. OOOE+OO 

1.91062E-01 

2 . 067E-01 

0. OOOE+OO 

2.70992E-01 

2 . 081E-01 

0. OOOE+OO 



3 . 59198E-01 
4.52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 
9.32816E-01 
9.72288E-01 
9 . 94700E-01 

Spatial point 0.20000 
Mu 

5 . 29953E-03 
2.77125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2.70992E-01 
3 . 59198E-01 
4.52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Spatial point 0.60000 
Mu 

5 . 29953E-03 
2 . 77 125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2.70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5 . 47506E-01 
6 . 40802E-01 
7 . 29008E-01 
8 . 08938E-01 


2 . 070E-01 

0 . 000E+00 

2 . 035E-01 

0 . 000E+00 

1 . 983E-01 

O.OOOE+OO 

1 . 924E-01 

0 . 000E+00 

1 . 864E-01 

0.000E+00 

1 . 808E-01 

0.000E+00 

1.760E-01 

0.000E+00 

1.722E-01 

0.000E+00 

1.696E-01 

0.000E+00 

1 . 681E-01 

0.000E+00 

Mu < 0 

Mu > 0 

2 . 187E-01 

2 . 081E-01 

2 . 142E-01 

2 . 089E-01 

2 . 056E-01 

2 . 106E-01 

1 . 937E-01 

2 . 123E-01 

1.799E-01 

2. 126E-01 

1 . 652E-01 

2 . 092E-01 

1 . 508E-01 

2 . 016E-01 

1.377E-01 

1 . 914E-01 

1 . 262E-01 

1 . 804E-01 

1 . 165E-01 

1 . 699E-01 

1 . 085E-01 

1 . 605E-01 

1 . 021E-01 

1 . 525E-01 

9.712E-02 

1.462E-01 

9 . 348E-02 

1.414E-01 

9 . 102E-02 

1.381E-01 

8 . 968E-02 

1 . 363E-01 

Mu < 0 

Mu > 0 

4 . 175E-02 

8.061E-02 

4 . 045E-02 

8 . 263E-02 

3 . 851E-02 

8 . 696E-02 

3.614E-02 

9.424E-02 

3 . 342E-02 

1 . 043E-01 

3.045E-02 

1 . 155E-01 

2.751E-02 

1.255E-01 

2.481E-02 

1.326E-01 

2 . 249E-02 

1 . 368E-01 

2 . 055E-02 

1.385E-01 

1 . 898E-02 

1 . 385E-01 

1 . 774E-02 

1 . 376E-01 



8.77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Spatial point 1.00000 
Mu 

5 . 29953E-03 
2.77125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4 . 52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Angular flux at the Direction 

Spatial point 0.00000 
Mu 

5 . 29953E-03 
2.77125E-02 
6 . 7 1844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4.52494E-01 
5.47506E-01 
6 . 40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9.94700E-01 

Spatial point 0.20000 


1 . 679E-02 

1 . 362E-01 

1 . 609E-02 

1.349E-01 

1.563E-02 

1.339E-01 

1 . 538E-02 

1 . 332E-01 

Mu < 0 

Mu > 0 

0 . 000E+00 

1 . 320E-02 

0 . 000E+00 

1 . 397E-02 

0 . 000E+00 

1 . 519E-02 

0 . 000E+00 

1.696E-02 

0.000E+00 

1 . 995E-02 

0 . 000E+00 

2 . 518E-02 

0 . 000E+00 

3 . 259E-02 

0 . 000E+00 

4 . 105E-02 

0 . 000E+00 

4.925E-02 

0 . 000E+00 

5 . 631E-02 

0 . OOOE+OO 

6. 192E-02 

0 . OOOE+OO 

6 . 610E-02 

0. 000E+00 

6 . 908E-02 

0. 000E+00 

7 . 108E-02 

0. 000E+00 

7 . 233E-02 

0. 000E+00 

7 . 297E-02 

Edit Points 

for group 3 

Mu <0 

Mu >0 

4 . 235E-02 

0. 000E+00 

4.407E-02 

0. 000E+00 

4 . 606E-02 

0. 000E+00 

4.823E-02 

0. 000E+00 

5 . 096E-02 

0. OOOE+OO 

5 . 385E-02 

0. OOOE+OO 

5 . 619E-02 

0. OOOE+OO 

5.773E-02 

0. OOOE+OO 

5 . 853E-02 

0. OOOE+OO 

5 . 878E-02 

0. OOOE+OO 

5 . 866E-02 

0. OOOE+OO 

5 . 834E-02 

0. OOOE+OO 

5 . 794E-02 

0. OOOE+OO 

5 . 756E-02 

0. OOOE+OO 

5 . 725E-02 

0. OOOE+OO 

5 . 707E-02 

0. OOOE+OO 
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Mu 

5 . 29953E-03 
2.77125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3.59198E-01 
4.52494E-01 
5.47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8 . 77702E-01 
9 . 32816E-01 
9.72288E-01 
9 . 94700E-01 

Spatial point 0.60000 
Mu 

5 . 29953E-03 
2 . 77125E-02 
6 . 7 1844E-02 
1 . 22298E-01 
1 . 91062E-01 
2 . 70992E-01 
3 . 59198E-01 
4. 52494E-01 
5 . 47506E-01 
6.40802E-01 
7 . 29008E-01 
8 . 08938E-01 
8.77702E-01 
9 . 32816E-01 
9 . 72288E-01 
9 . 94700E-01 

Spatial point 1.00000 
Mu 

5 . 29953E-03 
2 . 77 125E-02 
6.71844E-02 
1 . 22298E-01 
1 . 91062E-01 


Mu <0 

Mu > 0 

9 . 879E-02 

5.804E-02 

9.788E-02 

5.756E-02 

9 . 570E-02 

5 . 694E-02 

9 . 233E-02 

5 . 575E-02 

8 . 817E-02 

5 . 298E-02 

8.379E-02 

4 . 877E-02 

7 . 965E-02 

4.409E-02 

7 . 596E-02 

3 . 966E-02 

7 . 274E-02 

3 . 581E-02 

6 . 994E-02 

3 . 260E-02 

6 . 755E-02 

3 . 001E-02 

6 . 554E-02 

2 . 798E-02 

6.391E-02 

2 . 643E-02 

6 . 266E-02 

2 . 530E-02 

6.179E-02 

2 . 454E-02 

6. 131E-02 

2 . 413E-02 

Mu < 0 

Mu > 0 

5 . 391E-02 

4 . 129E-02 

5.371E-02 

4 . 145E-02 

5.324E-02 

4 . 217E-02 

5 . 249E-02 

4 . 367E-02 

5 . 147E-02 

4 . 601E-02 

5 . 021E-02 

4 . 875E-02 

4 . 873E-02 

5.107E-02 

4.707E-02 

5 . 247E-02 

4 . 533E-02 

5 . 294E-02 

4 . 360E-02 

5.271E-02 

4 . 199E-02 

5 . 204E-02 

4 . 057E-02 

5 . 120E-02 

3 . 939E-02 

5 . 035E-02 

3.847E-02 

4 . 962E-02 

3 . 783E-02 

4 . 907E-02 

3 . 747E-02 

4.875E-02 

Mu < 0 

Mu > 0 

0 . OOOE+OO 

1 . 815E-02 

0 . OOOE+OO 

1 . 885E-02 

0. 000E+00 

1 . 987E-02 

0. 000E+00 

2 . 115E-02 

0. 000E+00 

2 . 262E-02 
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2 . 70992E-01 

0 . OOOE+OO 

2.423E-02 

3 . 59198E-01 

0 . OOOE+OO 

2 . 592E-02 

4 . 52494E-01 

0. OOOE+OO 

2.759E-02 

5.47506E-01 

0. OOOE+OO 

2.916E-02 

6.40802E-01 

0. OOOE+OO 

3 . 057E-02 

7 . 29008E-01 

0. OOOE+OO 

3.174E-02 

8 . 08938E-01 

0. OOOE+OO 

3 . 266E-02 

8 . 77702E-01 

0. OOOE+OO 

3 . 335E-02 

9.32816E-01 

0. OOOE+OO 

3 . 383E-02 

9 . 72288E-01 

0. OOOE+OO 

3 . 414E-02 

9 . 94700E-01 

0. OOOE+OO 

3 . 430E-02 
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Unit 23 - group scalar fluxes in text format: 
Group Scalar Flux for group 1 


xb 

Flux 

0 . 00000E+00 

1 . 28432E+00 

2 . 00000E-02 

1 . 28890E+00 

4 . 00000E-02 

1 . 27715E+00 

6 . 00000E-02 

1 . 26014E+00 

8 . 00000E-02 

1 . 24036E+00 

1.00000E-01 

1 . 21900E+00 

1.20000E-01 

1 . 19682E+00 

1.40000E-01 

1 . 17436E+00 

1 . 60000E-01 

1 . 15213E+00 

1 . 80000E-01 

1 . 13071E+00 

2 . 00000E-01 

1 . 11225E+00 

2 . 20000E-01 

1 . 06075E+00 

2.40000E-01 

1 . 00266E+00 

2 . 60000E-01 

9.43901E-01 

2 . 80000E-01 

8 . 86170E-01 

3 . 00000E-01 

8 . 30328E-01 

3 . 20000E-01 

7 . 76878E-01 

3 . 40000E-01 

7 . 26119E-01 

3 . 60000E-01 

6 . 78219E-01 

3 . 80000E-01 

6 . 33262E-01 

4 . 00000E-01 

5 . 91277E-01 

4 . 20000E-01 

5 . 52258E-01 

4.40000E-01 

5 . 16184E-01 

4 . 60000E-01 

4 . 83027E-01 

4 . 80000E-01 

4 . 52768E-01 

5 . OOOOOE-Ol 

4 . 25412E-01 

5.20000E-01 

4 . 01012E-01 

5.40000E-01 

3.79713E-01 

5 . 60000E-01 

3 . 61861E-01 

5 . 80000E-01 

3 . 48338E-01 

6. OOOOOE-Ol 

3.43891E-01 

6 . 20000E-01 

3.43300E-01 

6.40000E-01 

3 . 37939E-01 

6 . 60000E-01 

3.30807E-01 

6 . 80000E-01 

3 . 22562E-01 

7. OOOOOE-Ol 

3 . 13515E-01 

7 . 20000E-01 

3 . 03852E-01 

7.40000E-01 

2 . 93694E-01 



7.60000E-01 

2 . 83123E-01 

7 . 80000E-01 

2.72198E-01 

8 . 00000E-01 

2.60959E-01 

8 . 20000E-01 

2.49433E-01 

8.40000E-01 

2 . 37633E-01 

8.60000E-01 

2 . 25560E-01 

8.80000E-01 

2 . 13198E-01 

9 . 00000E-01 

2 . 00514E-01 

9 . 20000E-01 

1 . 87448E-01 

9.40000E-01 

1 . 73885E-01 

9.60000E-01 

1 . 59609E-01 

9 . 80000E-01 

1.44114E-01 

1 . OOOOOE+OO 

1 . 24690E-01 

Group Scalar Flux for group 

xb 

Flux 

0. 00000E+00 

1 . 94270E-01 

2 . 00000E-02 

2 . 51740E-01 

4 . 00000E-02 

2 . 83754E-01 

6 . 00000E-02 

3 . 05655E-01 

8 . 00000E-02 

3 . 20727E-01 

1 . 00000E-01 

3 . 30516E-01 

1 . 20000E-01 

3 . 35915E-01 

1.40000E-01 

3 . 37475E-01 

1 . 60000E-01 

3 . 35525E-01 

1 . 80000E-01 

3.30200E-01 

2 . 00000E-01 

3 . 21957E-01 

2 . 20000E-01 

3 . 17179E-01 

2.40000E-01 

3 . 10880E-01 

2 . 60000E-01 

3 . 03654E-01 

2 . 80000E-01 

2 . 95773E-01 

3 . 00000E-01 

2 . 87427E-01 

3 . 20000E-01 

2 . 78751E-01 

3.40000E-01 

2 . 69845E-01 

3 . 60000E-01 

2 . 60789E-01 

3 . 80000E-01 

2.51648E-01 

4 . 00000E-01 

2.42470E-01 

4 . 20000E-01 

2 . 33295E-01 

4 . 40000E-01 

2 . 24150E-01 

4 . 60000E-01 

2 . 15053E-01 

4 . 80000E-01 

2 . 06008E-01 

5 . 00000E-01 

1 . 97009E-01 

5 . 20000E-01 

1 . 88030E-01 
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5.40000E-01 1.79015E-01 

5 . 60000E-01 1.69851E-01 

5.80000E-01 1.60282E-01 

6 . 00000E-01 1.49053E-01 

6 . 20000E-01 1 . 38195E-01 

6.40000E-01 1 . 29726E-01 

6 . 60000E-01 1 . 22260E-01 

6 . 80000E-01 1 . 15486E-01 

7 . 00000E-01 1 . 09246E-01 

7 . 20000E-01 1 . 03439E-01 

7 . 40000E-01 9 . 79969E-02 

7 . 60000E-01 9 . 28664E-02 

7 . 80000E-01 8 . 80063E-02 

8 . 00000E-01 8 . 33824E-02 

8.20000E-01 7 . 89658E-02 

8 . 40000E-01 7.47310E-02 

8 . 60000E-01 7 . 06547E-02 

8 . 80000E-01 6 . 67145E-02 

9 . 00000E-01 6 . 28879E-02 

9 . 20000E-01 5.91494E-02 

9 . 40000E-01 5 . 54677E-02 

9.60000E-01 5 . 17967E-02 

9 . 80000E-01 4 . 80459E-02 

1 . OOOOOE+OO 4 . 37958E-02 
Group Scalar Flux for group 3 
xb Flux 

0. 00000E+00 5 . 52502E-02 

2 . 00000E-02 6 . 58401E-02 

4 . 00000E-02 7 . 31773E-02 

6 . 00000E-02 7 . 93078E-02 

8 . 00000E-02 8.46920E-02 

1.00000E-01 8 . 95801E-02 

1 . 20000E-01 9 . 41610E-02 

1 . 40000E-01 9 . 86202E-02 

1 . 60000E-01 1 . 03195E-01 

1.80000E-01 1 . 08321E-01 

2 . 00000E-01 1 . 16030E-01 

2 . 20000E-01 1 . 2267 IE-01 

2 . 40000E-01 1 . 25235E-01 

2 . 60000E-01 1 . 26044E-01 

2 . 80000E-01 1 . 25716E-01 

3 . 00000E-01 1 . 24596E-01 



3 . 20000E-01 
3 . 40000E-01 
3.60000E-01 
3 . 80000E-01 
4.00000E-01 
4.20000E-01 
4.40000E-01 
4 . 60000E-01 
4.80000E-01 
5 . 00000E-01 
5 . 20000E-01 
5 . 40000E-01 
5 . 60000E-01 
5 . 80000E-01 
6 . 00000E-01 
6 . 20000E-01 
6 . 40000E-01 
6 . 60000E-01 
6 . 80000E-01 
7 . 00000E-01 
7 . 20000E-01 
7 . 40000E-01 
7 . 60000E-01 
7 . 80000E-01 
8 . 00000E-01 
8 . 20000E-01 
8.40000E-01 
8 . 60000E-01 
8 . 80000E-01 
9 . 00000E-01 
9.20000E-01 
9.40000E-01 
9 . 60000E-01 
9.80000E-01 
1 . OOOOOE+OO 


1 . 22911E-01 
1.20823E-01 
1 . 18458E-01 
1 . 15917E-01 
1 . 13282E-01 
1 . 10622E-01 
1.07998E-01 
1 . 05463E-01 
1 . 03065E-01 
1 . 00853E-01 
9 . 88782E-02 
9 . 71995E-02 
9 . 58986E-02 
9 . 51256E-02 
9 . 55665E-02 
9 . 58034E-02 
9.44263E-02 
9 . 24749E-02 
9 . 01703E-02 
8 . 76107E-02 
8.48520E-02 
8 . 19290E-02 
7 . 88645E-02 
7 . 56735E-02 
7 . 23646E-02 
6 . 89415E-02 
6 . 54029E-02 
6 . 17417E-02 
5 . 79443E-02 
5 . 39875E-02 
4 . 98336E-02 
4 . 54193E-02 
4 . 0627 IE-02 
3 . 51899E-02 
2.77138E-02 
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Unit 24 - total scalar flux in text format: 

Total Scalar Flux 

x Flux 

0 . 00000E+00 1 . 53384E+00 

2.00000E-02 1.60648E+00 

4.00000E-02 1 . 63408E+00 

6.00000E-02 1 . 64511E+00 

8 . 00000E-02 1 . 64578E+00 

1 . 00000E-01 1 . 63910E+00 

1 . 20000E-01 1 . 62689E+00 

1.40000E-01 1 . 61046E+00 

1 . 60000E-01 1 . 59085E+00 

1 . 80000E-01 1 . 56923E+00 

2 . 00000E-01 1 . 55024E+00 

2 . 20000E-01 1 . 50061E+00 

2.40000E-01 1 . 43877E+00 

2 . 60000E-01 1 . 37360E+00 

2 . 80000E-01 1 . 30766E+00 

3 . 00000E-01 1 . 24235E+00 

3 . 20000E-01 1 . 17854E+00 

3 . 40000E-01 1 . 11679E+00 

3 . 60000E-01 1.05747E+00 

3 . 80000E-01 1 . 00083E+00 

4 . 00000E-01 9.47028E-01 

4 . 20000E-01 8 . 96176E-01 

4.40000E-01 8.48333E-01 

4. 60000E-01 8.03543E-01 

4.80000E-01 7 . 61842E-01 

5 . 00000E-01 7 . 23275E-01 

5 . 20000E-01 6 . 87920E-01 

5.40000E-01 6 . 55928E-01 

5 . 60000E-01 6 . 27611E-01 

5.80000E-01 6 . 03746E-01 

6 . 00000E-01 5 . 88511E-01 

6 . 20000E-01 5 . 77298E-01 

6.40000E-01 5 . 62092E-01 

6 . 60000E-01 5.45542E-01 

6 . 80000E-01 5 . 28218E-01 

7 . OOOOOE-Ol 5 . 10371E-01 
7 . 20000E-01 4 . 92144E-01 

7.40000E-01 4 . 73620E-01 



7 . 60000E-01 
7 . 80000E-01 
8.00000E-01 
8 . 20000E-01 
8.40000E-01 
8 . 60000E-01 
8 . 80000E-01 
9 . 00000E-01 
9 . 20000E-01 
9 . 40000E-01 
9 . 60000E-01 
9 . 80000E-01 
1 . OOOOOE+OO 


4 . 54854E-01 
4.35878E-01 
4. 16706E-01 
3 . 97341E-01 
3.77767E-01 
3.57956E-01 
3 . 37857E-01 
3 . 17390E-01 
2.96431E-01 
2.74772E-01 
2 . 52033E-01 
2.27350E-01 
1 . 96200E-01 
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